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In  this  paper,  we  construct  very  efficient  high-order  schemes  for  general  time-dependent  advection- 
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derivatives  of  the  source  terms.  From  the  latter  technique,  we  construct  spatially  third-,  fourth-,  and 
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Published  by  Elsevier  Ltd. 


1.  Introduction 

In  this  paper,  we  construct  very  efficient  high-order  schemes  for 
general  time-dependent  advection-diffusion  problems,  based  on 
the  first-order  hyperbolic  system  method  [1,2],  In  this  method, 
the  diffusion  term  is  reformulated  as  a  hyperbolic  system,  leading 
to  the  unification  of  advection  and  diffusion  as  a  single  hyperbolic 
system  [2],  The  drastic  change  in  the  type  of  equations,  parabolic  to 
hyperbolic,  brings  several  dramatic  improvements  in  the  construc¬ 
tion  of  numerical  schemes:  hyperbolic  schemes  for  diffusion,  the 
same  order  of  accuracy  for  the  solution  and  the  gradients, 
orders-of-magnitude  convergence  acceleration,  etc.,  which  have 
been  demonstrated  for  steady  diffusion  and  viscous  problems  in 
Refs.  [1-5]  and  unsteady  advection-diffusion  problems  in  Ref. 
[6],  It  is  based  on  the  reformulation  of  the  governing  equations, 
and  therefore  applicable  to  any  discretization  method.  In  this  work, 
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we  consider  a  Residual-Distribution  (RD)  method  [7],  which  has 
been  well  developed  for  hyperbolic  systems  and  has  a  superior  fea¬ 
ture  of  achieving  second-order  accuracy  in  a  compact  stencil. 

In  the  previous  work  [6],  we  extended  the  hyperbolic  method, 
for  the  first  time,  to  time-accurate  computations  by  an  implicit 
time-integration  method  based  on  the  second-order  backward 
difference  formula.  The  resulting  scheme  was  applied  to  various 
time-dependent  problems,  demonstrating  second-order  accuracy 
for  the  solution  and  the  gradient  achieved  at  all  interior  and 
boundary  nodes  in  uniform  and  nonuniform  grids  at  every  physical 
time  step,  and  rapid  convergence  for  solving  implicit-residual 
equations  by  Newton’s  method  (i.e.,  less  than  5  iterations  per 
physical  time  step),  which  is  possible  by  the  compactness  of  the 
RD  schemes.  As  a  consequence  of  the  first-order  re-formulation 
of  the  equation,  the  number  of  linear  relaxations  performed  at 
every  Newton  iteration  was  shown  to  increase  only  linearly  with 
the  grid  size,  not  quadratically  as  typical  for  diffusion  problems. 
The  efficiency  of  the  developed  second-order  schemes  was  demon¬ 
strated  for  linear  and  nonlinear  advection-diffusion  problems  on 
highly  refined  grids,  up  to  30,000  nodes. 
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In  this  paper,  we  propose  a  very  simple  extension  of  the  second- 
order  schemes  to  higher-order.  We  show  that  high-order  spatial 
accuracy  can  be  achieved  simply  by  modifying  the  source  term 
discretization.  There  are  two  approaches  to  the  source  term 
discretization:  (1)  reformulation  of  the  source  terms  with  their 
divergence  forms  and  (2)  correction  to  the  trapezoidal  rule  for 
the  source  term  discretization.  The  former  technique  is  based  on 
the  divergence  formulation  of  source  terms  proposed  in  Ref.  [8]: 
write  the  source  term  in  the  divergence  form  and  discretize  it  in 
the  same  way  as  the  flux  divergence  term.  The  latter  is  based  on 
a  high-order  correction  to  the  trapezoidal  rule,  and  thus  called  here 
the  generalized  trapezoidal  rule.  In  either  case,  high-order  accu¬ 
racy  is  achieved  by  making  low-order  truncation  error  terms  pro¬ 
portional  to  the  residual,  which  thus  vanish  in  the  steady  state 
and  yield  high-order  accuracy.  We  solve  the  resulting  implicit- 
residual  equations  by  an  implicit  solver  based  on  the  second-order 
Jacobian  matrix  developed  in  the  previous  work  [6],  As  we  will 
show,  the  implicit  solver  is  as  powerful  as  Newton’s  method;  e.g., 
eight  orders  of  magnitude  reduction  can  be  achieved  in  10  itera¬ 
tions.  To  enable  time-accurate  computations,  we  employ  high- 
order  versions  of  the  backward  difference  formulas  (BDF),  which 
are  incorporated  as  source  terms,  and  solve  the  implicit-residual 
equations  by  the  implicit  solver  over  each  physical  time  step.  In 
this  manner,  the  steady  state  is  made  equivalent  to  the  next 
physical  time  with  all  the  benefits  of  the  hyperbolic  method 
retained.  We  note  that  the  choice  of  the  implicit  time  stepping 
method  is  independent  of  the  developed  high-order  RD  schemes, 
and  thus  other  methods  such  as  implicit  Runge-Kutta  methods 
or  space-time  methods  can  also  be  employed. 

The  high-order  RD  schemes  developed  in  this  work  are  signifi¬ 
cantly  different  from  other  high-order  RD  schemes  in  that  our 
schemes  are  based  on  the  first-order  hyperbolic  system  formula¬ 
tion  of  the  advection-diffusion  equation  2].  In  this  approach,  the 
loss  of  high-order  accuracy  in  the  intermediate  Reynolds  number, 
as  discussed  in  Refs.  [9-11],  cannot  occur  because  the  advective 
and  diffusive  terms  are  fully  integrated  into  a  single  hyperbolic 
system.  If  the  original  advection-diffusion  equation  is  discretized, 
a  high-order  RD  scheme  needs  to  be  developed  for  the  diffusion 
term  (i.e.,  second  derivative)  and  then  carefully  combined  with 
an  advection  scheme,  e.g.  by  using  a  blending  parameter  as 
described  in  Ref.  [10],  to  avoid  the  loss  of  accuracy.  Furthermore, 
while  high-order  RD  schemes  based  on  high-order  elements 
require  extra  degrees  of  freedom  for  each  variable,  our  schemes 
are  based  on  linear  elements  for  any  order  of  accuracy  but  require 
extra  gradient  variable  to  be  added  to  the  solution  vector.  Note  that 
the  number  of  extra  variables  in  the  high-order  elements  increase 
for  higher-order  accuracy,  but  the  number  of  extra  variables 
required  in  our  approach  is  fixed  and  independent  of  the  order  of 
accuracy.  Our  approach  is  somewhat  similar  to  those  in  Refs. 
[12-14],  but  again  is  significantly  different  by  the  use  of  first-order 
hyperbolic  system  formulation  of  the  advection-diffusion  equation 
and  by  the  source  term  discretization  techniques.  It  is  emphasized 
that  our  schemes  require  only  the  first  and  second  derivatives  of 
the  source  term,  or  in  some  cases  the  first  derivatives  only;  they 
do  not  require  the  gradient  computation  for  the  solution  variables. 

The  third-order  schemes  developed  in  this  paper  are  similar  to 
the  third-order  finite-volume  scheme  of  Katz  and  Sankaran  [15,16] 
in  that  the  second-order  truncation  error  is  eliminated  by  making 
it  proportional  to  the  residual  and  the  upgrade  is  achieved  by  sec¬ 
ond-order  accurate  gradients.  However,  as  we  demonstrate  in  this 
paper,  the  proposed  high-order  RD  schemes  have  several  superior 
features:  (1)  implicit  solver  can  be  constructed  by  the  Jacobian 
derived  from  a  compact  second-order  RD  scheme,  (2)  gradient 
computations  are  required  for  the  source  terms  only,  and  not  for 
the  solution,]  3)  stiffness  due  to  the  second  derivative  of  the  diffu¬ 
sion  term  is  completely  eliminated,  (4)  higher-order  schemes  can 


be  constructed  beyond  third-order  (in  extending  it  to  multi¬ 
dimensions),  and  (5)  the  same  order  of  accuracy  is  achieved  for 
the  gradients,  as  well.  In  particular,  the  fourth-order  scheme 
constructed  in  Section  5  is  remarkably  more  efficient  because  it 
does  not  require  second  derivatives  of  the  source  term,  which  are 
required  in  the  schemes  described  in  Refs.  [15,16], 

In  this  paper,  we  focus  on  one-dimensional  linear  and  nonlinear 
advection-diffusion  problems.  It  certainly  serves  as  a  basis  for  the 
development  of  high-order  multi-dimensional  RD  schemes  for 
more  complex  equations.  Yet,  more  importantly,  the  one-dimen¬ 
sional  high-order  schemes  developed  in  this  paper  could  poten¬ 
tially  bring  significant  improvements  to  practical  problems  such 
as  material  thermal  response  calculations  of  thermal  protection 
systems  of  atmospheric  entry  vehicles  [17-19],  and  the  experi¬ 
mental  aeroheating  data  reduction  [20,21],  which  are  based  on 
one-dimensional  analyses  and  routinely  used  in  industries  (e.g. 
NASA).  The  extension  to  higher  dimensions  is  beyond  the  scope 
of  the  paper;  it  will  be  addressed  in  a  subsequent  paper. 

The  paper  is  organized  as  follows.  In  the  next  section,  the  time- 
dependent  hyperbolic  advection-diffusion  system  is  described.  In 
Section  3,  a  compact  second-order  residual-distribution  scheme, 
a  steady  solver,  and  the  second-order  discretization  are  discussed. 
In  Section  4,  the  third-  and  fourth-order  RD  schemes  with  source 
term  reformulation  are  proposed.  In  Section  5,  the  third-,  fourth, 
and  sixth-order  RD  schemes  with  source  term  discretization  are 
developed  and  proposed.  Numerical  results  are  then  presented  in 
Section  6.  Finally,  Section  7  concludes  the  study  with  remarks. 

2.  Time-dependent  hyperbolic  advection-diffusion  system 

We  start  with  a  linear  advection-diffusion  equation  to  simplify 
the  discussion.  We  will  extend  the  discussion  later  to  a  more  gen¬ 
eral  nonlinear  advection-diffusion  equation. 

Consider  the  one-dimensional  (1-D)  time-dependent  advec¬ 
tion-diffusion  equation: 

dtu  +  adxu  =  vdxxii  +  S(x),  (1) 

where  a  and  v  are  both  taken  to  be  positive  constant,  and  S  is  the 
source  term.  We  will  follow  the  procedure  we  described  in  Ref. 
[6]  and  rewrite  the  above  equation  as  a  first-order  hyperbolic 
advection-diffusion  system: 

dTu  = -adxU  +  Vdxp-^u  +  Six),  (2) 

d Tp  =  (dxu  -  p)/Tr,  (3) 

where  the  relaxation  time,  Tr  >  0,  is  arbitrary  and  defined  as 
described  in  Ref.  [6],  and  S  includes  any  existing  source  terms  from 
the  advection-diffusion  problem,  S,  as  well  as  any  additional  terms 
that  arise  from  the  implicit  time-stepping  scheme,  At  is  the  physical 
time  steps,  and  z  is  the  pseudo  time  step.  Note  that  the  dtp  is  taken 
as  pseudo  time  derivative,  dzp. 

The  variable  a  depends  on  the  order  of  the  Backward-Differenc- 
ing-Formula  (BDF):  1  for  the  lst-order  (BDF1),  3/2  for  the  second- 
order  (BDF2),  11/6  for  the  third-order  (BDF3),  25/12  for  the  fourth- 
order,  and  147/60  for  the  sixth-order  time  discretizations  (see 
Table  1).  The  remaining  terms  in  the  BDF  are  stored  in  the  source 
term  function  S(x).  It  is  well  known  that  the  BDF2  is  A-stable  and 
higher-order  BDFs  are  not.  Therefore,  the  second-order  scheme  is 
unconditionally  stable,  but  higher-order  BDFs  are  conditionally 
stable.  Consequently,  the  stability  of  the  higher-order  schemes 
depends  on  the  spatial  discretization.  Estimates  for  the  maxi¬ 
mum-allowable  CFL  numbers  are  given  in  Appendix  A  for  a  set  of 
representative  high-order  schemes  developed  in  this  paper. 

Towards  the  pseudo  steady  state,  the  variable  p  approaches  the 
solution  gradient  and  hence  the  above  equation  becomes  identical 
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Table  1 

BDF  coefficients. 


BDF  order 

un+l 

Un 

U"-1 

u"-2 

u„-3 

u"-4 

un-5 

1 

1 

-1 

2 

3 

-2 

1 

3 

11 

-3 

3 

1 

4 

25 

12 

-4 

3 

4 

3 

1 

4 

5 

137 

5 

5 

10 

5 

1 

60 

3 

4 

5 

6 

147 

-6 

15 

20 

15 

6 

1 

60 

T 

T 

T 

3 

6 

to  the  original  advection-diffusion  equation  with  the  time  deriva¬ 
tive  discretized  by  the  BDF,  i.e.,  a  consistent  discretization  of  the 
original  equation.  The  system  reduces,  of  course,  to  the  original 
steady  equation  also  in  the  physical  steady  state  when  dcu  =  0. 
Second-order  accurate  unsteady  computations  have  been 
demonstrated  based  on  the  above  formulation  in  Ref.  [6]. 

In  the  vector  form,  our  time-dependent  first-order  advection- 
diffusion  system  can  be  written  as 


<9U  SU 
<9t  +  dx 


=  S, 


(4) 


where 


U  = 


'-au/At  +  S(x)' 

.  ~P/Tr  . 


(5) 


For  any  positive  Tr,  the  Jacobian  has  the  following  two  real 
eigenvales: 


Ai 


1-1  1  + 


4v 


a 

2 


1  + 


(6) 


with  linearly  independent  eigenvectors  [6],  and  therefore  the  above 
system  is  hyperbolic  in  the  pseudo  time.  The  hyperbolicity  serves 
mainly  as  a  guide  for  discretization:  various  discretization  tech¬ 
niques  are  available  for  hyperbolic  systems,  e.g.,  upwinding.  In 
addition  to  the  convenience  in  discretization,  the  major  benefits 
are:  (1)  the  hyperbolic  discretization  results  in  a  strong  coupling 
between  the  two  variables  that  results  in  the  equal  order  of  accu¬ 
racy  for  both  u  and  p  =  dxu  on  arbitrary  grids  and  (2)  0(1 /h)  accel¬ 
eration  is  achieved  in  iterative  convergence  for  the  linearized 
residual  equation  in  implicit  solvers  due  to  the  0(1  /h)  condition 
number  of  the  Jacobian,  which  is  0(h)  smaller  than  that  of  the  Jaco¬ 
bian  derived  from  a  conventional  diffusion  scheme.  For  explicit 
schemes,  the  0(1  /h)  acceleration  is  achieved  by  the  0(h)  time  step 
typical  in  methods  for  hyperbolic  systems,  which  is  0(  1  /h)  times 
larger  than  a  typical  time  step  of  0(h2)  for  diffusion.  The  0(1 /h) 
acceleration  in  convergence  has  been  demonstrated  over  traditional 
methods  for  the  diffusion  equation  [1,4,  for  the  advection-diffusion 
equation  [2,5],  for  the  compressible  Navier-Stokes  equations  [3], 
and  most  recently  for  time-dependent  linear  and  nonlinear  advec¬ 
tion-diffusion  problems  using  RD  technique  [6[. 

In  the  next  sections,  we  first  briefly  describe  the  second-order 
time-dependent  discretization  process  using  the  RD  scheme.  Fur¬ 
ther  details  are  given  in  Ref.  [6[.  We  then  extend  the  order  of  accu¬ 
racy  of  the  scheme  to  third-order  and  fourth-order  RD  hyperbolic 
advection-diffusion  schemes  with  reformulation  of  the  hyperbolic 
system,  and  finally  to  higher  order  (fourth-  and  sixth-order)  RD 
hyperbolic  advection-diffusion  schemes  with  introduction  of 
new  source  term  discretization  technique. 


3.  Second-order  RD  hyperbolic  advection-diffusion 

The  RD  method  requires  (1)  evaluation  of  the  cell  (or  element) 
residuals  and  (2)  the  distribution  of  the  residuals  to  the  nodes 


bounding  the  cell.  Consider  a  one-dimensional  domain  discretized 
with  N  nodes  that  are  distributed  arbitrary  over  the  domain  of 
interest  with  the  solution,  u,  and  the  solution  gradient,  p  =  dxu, 
data  stored  at  the  nodes  denoted  by  Xj,  i  =  1,2,3, ...  ,N  (Fig.  1). 
We  define  the  cell-residual  <bc  by  integrating  the  spatial  part  of 
Eq.  (4)  over  the  cell,  C,  defined  by  the  nodes  i  and  i  +  1 : 


rx i+i 

:  /  (-AUx  +  S)dx, 

JXj 

-a(ui+ 1  -  Ui)  +  v(pi+1  -  pi, 

Tr  [Ui+1  -  “01 


(7) 

(8) 


The  first  term  of  Eq.  (8)  is  the  result  of  the  exact  integration  of 
AU„  over  the  cell  C.  The  integration  of  the  second  term  is  not 
exact  and  therefore  is  the  source  of  the  overall  discretization  error. 
We  will  further  discuss  this  discretization  error  in  the  following 
sections. 

We  construct  the  spatial  discretization  of  Eq.  (4)  by  distributing 
the  cell  residuals  to  the  nodes: 


^1  =  1  (B+<Kl  +  fT<f>R),  (9) 

where  <t>'  and  <PR  denote  the  cell-residuals  over  the  left  and  right 
cells  of  the  node  i,  respectively,  and  hi  is  the  dual  volume  (see 
Fig.  1 )  defined  by 

^  =  hL  +  hR .  hL  =  Xj  -  x,  ! ,  hR  =  xi+i  -  Xj.  (10) 

Note  that  the  pseudo  time  derivative  on  the  left  hand  side  is 
retained  here  just  for  the  sake  of  illustration.  It  will  be  ignored  in 
the  implicit  formulation  that  follows.  Our  aim  is  to  directly  solve 
it  for  the  pseudo  steady  state,  which  corresponds  to  the  next  phys¬ 
ical  time.  In  this  sense,  our  method  is  not  a  dual-time  stepping 
method.  The  distribution  matrices  £T  and  B+  (Fig.  2)  typically  do 
not  affect  the  order  of  accuracy  of  the  discretization  scheme;  it  is 
determined  by  the  cell-residuals  [22,23],  and  therefore  our  focus 
will  be  on  the  residual  evaluation.  (Refer  to  Ref.  [6]  for  formulation 
of  the  distribution  matrices  B~  and  B+). 


3.1.  Second-order  discretization 


We  may  use  a  simple  trapezoidal  rule  for  the  integration  of  the 
source  term  S  over  the  cell  C  and  arrive  at  the  following  cell 
residuals: 


<PL  = 


-a(ui+ 1  -  +  v(pM  -  p^  -  (xi+i  -  Xi)ft(ui+ 1  +  u^/ 2 

(“i+i  ~  ^  -  (xi+i  -  Xi)(pM  +Pi)/2)/Tr 
(xi+i  -Xi)(sj+i  +s/)/2 

0 


fe+i 


(11) 


where  k  and  n  are  the  Newton  iteration  counter  (as  described 
below)  and  the  physical  time  index,  respectively.  Note  that  the 


i-1 


K 


hR 


-o- 

i+1 


Fig.  1.  Schematic  of  grid  spacing  for  a  1-D  grid. 


S+0;_!  £“<£,  B+<D,  B«D,+1 

-o - - — ^-o^— - o- 

i-1  i  i+1  i+2 

Fig.  2.  Residual  distribution  to  the  nodes. 
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second  term  of  Eq.  (11),  which  is  computed  at  the  two  previous 
physical  time  steps,  is  constant  during  the  Newton  iteration,  and 
thus  will  not  contribute  to  the  Jacobian. 

The  implicit  formulation  is  defined  by 

Uk+1  =  Uk  +  AU\  (12) 


where  U  =  (ui,p1,u2,p2 _ ,uN,pN)  and  k  is  the  iteration  counter. 

The  correction  AUk  =  Uk+1  -  Uk  is  determined  as  the  solution  to 
the  linear  system: 


<9Res 

OU 


AUk 


=  -Res'1, 


(13) 


where  Resk  is  the  right  hand  side  of  Eq.  (9),  which  is  the  unsteady 
residual  vector  evaluated  by  Uk.  Note  that  the  pseudo  time 
derivative  has  been  ignored.  The  Jacobian  matrix  is  sparse,  having 
three  2x2  blocks  in  each  row  for  all  interior  nodes  and  two  blocks 
for  boundary  nodes.  The  ith  pair  of  rows  of  the  linear  system  is 
given  by 


Ji-i&UU  +J,  AU?  +JMAUk+i  =  -1  (B+<PL  +  B-<PR)\ 

(14) 

where  Al/f_1  =  (Auf , , Apf ,),  AUk  =  (A uf,  Apf),  Alfk+1 

=  (Au?+i,Apf+1), 

1  d(B+<PL) 

J hi  dUt ’ 

(15) 

1  fd(B+<PL)  (dB-<PR)\ 

Jl  h,  ^  dUi  +  dUt  )’ 

(16) 

1  d(B  <t>R) 

J,+1  ^  dUM  ■ 

(17) 

We  note  that  the  derivative  of  the  distribution  matrices,  IT  and  B+, 
are  zero  for  linear  advection-diffusion  problems  and  for  nonlinear 
problems  with  constant  a/v  values,  which,  for  simplicity,  are  the 
only  nonlinear  system  considered  in  this  paper.  However,  the 
proposed  schemes  are  applicable  to  general  nonlinear  advection- 
diffusion  problems. 

The  residual  Jacobians  needed  in  the  Newton  solver  are  exactly 
in  the  same  form  as  the  above  equations,  but  the  derivatives  of  the 
cell-residuals  now  include  the  contribution  from  the  physical  time 
derivative: 


d<PR 

dUt 


d  <PL 
d  U, 


a-h*m 

1 

-V  ‘ 
tlR 

(18) 

Tr 

2Tr . 

~a~hLj£-t 

V 

_  A 

(19) 

2  Tr 


At  each  physical  time  level  n,  we  solve  the  pseudo  steady  problem 
by  Newton’s  method  with  the  current  solution  at  n  as  the  initial 
solution.  This  results  in  second-order  discretization,  which  was 
demonstrated  with  several  examples  (linear  and  nonlinear)  in  Ref. 

[6], 

We  now  focus  on  the  truncation  error  of  the  source  term  dis¬ 
cretization  as  it  is  needed  for  our  discussion  in  the  next  sections 
where  we  show  two  different  approaches  to  obtain  higher-order 
RD  schemes.  We  demonstrate  the  truncation  error  of  the  system 
by  considering  the  residual  of  the  second  equation  in  our  hyper¬ 
bolic  advection-diffusion  system: 

<2>p  =  (Ui+i  -  u,  -  h(pi+1  +  p{)/2)/Tr,  (20) 

where  h  is  the  width  of  the  cell  C  (i.e.  xi+i  -  x,).  We  now  expand  the 
$p  around  the  node  i  to  obtain  the  truncation  error  (T.E.)  for  the  sec¬ 
ond  equation,  dzp,  which  after  some  algebra  and  rearrangements 
becomes: 


T.E.(drp)  = 


(• dxUi  -  Pi)  +  ^  (dxxUi  -  dxPi) 


+  -Q  {dxxxlli  ~  ^dxxPi)  +  0(h  ) 


Tr  =  0(h2), 


(21) 


where  the  first  two  terms  will  simply  vanish  as  they  satisfy  the  ori¬ 
ginal  equation  in  the  pseudo  steady  state.  The  non-vanishing  last 
term  makes  the  current  discretization  second-order.  The  same 
results  are  obtained  with  the  first  equation  and  therefore  is  not 
repeated  here.  The  above  analysis  shows  that  the  error  arises  from 
the  source  term  discretization,  not  from  the  flux  balance  term, 
which  is  exact  in  one  dimension.  It  implies  that  we  can  achieve 
high-order  accuracy  only  by  improving  the  source  term  discretiza¬ 
tion.  In  fact,  in  Ref.  [14],  a  fourth-order  finite-difference  scheme  is 
constructed  in  this  manner,  i.e.,  by  a  high-order  source  term  discret¬ 
ization  with  explicit  pseudo-time  stepping  targeted  for  steady 
problems. 

Note  that  the  above  truncation  error  analysis  is  based  on  the 
cell-residual  expanded  at  a  node.  Hence,  the  nodal  residual  has 
the  same  order  as  the  cell-residuals  because  it  is  constructed  as  a 
weighted  sum  of  the  cell-residuals  as  shown  in  Eq.  (9).  Note  also 
that  the  truncation  error  analysis  based  on  the  cell-residual  is  local 
to  the  cell  and  thus  valid  for  any  size  of  the  cell,  i.e.,  valid  for  uni¬ 
form  and  nonuniform  grids,  in  the  rest  of  the  paper,  the  truncation 
error  analysis  is  all  based  on  the  cell-residual. 

in  the  next  sections,  we  show  two  different  techniques  that  will 
lead  to  higher-order  schemes.  In  particular,  we  discuss  on  ( 1 )  refor¬ 
mulation  of  the  source  term  with  its  divergence  form  and  (2)  new 
source  term  discretization  technique  which  acts  as  a  correction  to 
the  trapezoidal  rule. 


4.  High-order  RD  scheme  with  source  term  reformulation 
(third-  and  fourth-order) 

Consider  the  source  term,  S,  as  the  divergence  of  a  function 
fs  -fx  =  S.  We  now  replace  the  source  terms  with  their  divergence 
forms,  and  rewrite  our  first-order  hyperbolic  advection-diffusion 
equation  as 


dzu  =  - adxu  +  vdxp  +  8*% ,  (22) 

dxp  =  (dxu  -  dJl)/Tr,  (23) 

where  /J  and  fjj  are  the  divergence  formulation  of  the  source  terms 
in  the  u  and  the  p  equations,  respectively.  With  the  above  reformu¬ 
lation  of  the  advection-diffusion  equation,  which  will  be  discussed 
in  more  details,  the  residual  evaluations  of  the  system  becomes 
exact  with  no  special  discretization  scheme: 


<PC=  r+,(-AUx  +  f^)dx,  (24) 

JXj 

(  -a(uM  -  uf)  +  v(pM  -  p^  +  (ff)m  -  (fus)i 
=  <  r  i  ■  (25) 

1  h  ["hi  -  -  C fp)i+1  +  (£)f)] 

We  remark  that  even  though  the  residual  evaluation  of  our  refor¬ 
mulated  advection-diffusion  system  is  exact,  the  overall  accuracy 
of  the  scheme  depends  on  the  accuracy  of  the  divergence  formula¬ 
tion  of  the  source  terms  and  how  it  is  formulated. 


4.1.  Third-order  scheme  with  divergence  form 

Following  the  divergence  formulation  introduced  in  [8],  we 
write  the  source  flux  in  a  more  general  form: 
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m> 3  / _ “|  \n-l 

-'S’  (26) 

n=l 

=  (X  -  Xi)S  -  ^  (x  -  Xi)2dxS  +  -L  (x  -  XifdxxS 

-^(x-xi)4dxxxS  +  ---,  (27) 

where  the  source  term  S  and  its  derivatives  Sx ,  S**,  etc.  are  not 
constants  but  functions  that  vary  in  space.  We  remark  that  the 
source  flux,  fs,  is  not  necessarily  a  polynomial  of  order  m  because 
it  depends  on  the  derivatives  of  the  source  term  S.  In  this  paper, 
we  do  not  consider  high-order  schemes  that  require  evaluation  of 
third  or  higher  derivatives.  The  divergence  of  the  source  flux,  djs, 
for  the  third-order  accuracy  (i.e.  m  =  3)  is 

t  h3  3 

djs  =  S  +  -dxxxS  =  S  +  0(h3),  (28) 

which  is  identical  to  the  original  equation  up  to  the  third-order.  We 
remark  that  this  error  does  not  necessarily  limit  the  maximum 
order  attained  by  numerical  schemes.  In  the  next  section,  we  will 
show  that  a  fourth-order  scheme  can  be  constructed  by  a  slightly 
modified  divergence  form  with  m  =  3. 

Discretization  of  the  reformulated  hyperbolic  system  leads  to 
the  cell  residual,  for  example,  for  the  second  equation  as 

<Pcp  =  [ui+1  -  ut  -  (fp\+1  +  (£).] /tTi  (29) 

where  special  care  must  be  taken  when  we  evaluate  /js  and  fZl , 
(because  of  the  presence  ofx,  in  the  divergence  formulation)  as  they 
depend  on  whether  the  cell  residual  is  being  distributed  to  the  left 
or  to  the  right  node  of  the  cell  (Fig.  3): 


=  [ui+1  -  u,-  -  C^)i+1  + 

h2  h3 

tf  |  I  ““  t/j  —  IIrSi  |  1  -p  ()ySt  I  I  ^  dxxSj  |  1 


TV, 


(30) 


0 


C, 

P 


[uM-u,-(f3)M  +  (f^/rr, 

Ui+i  —  U;  +  hftSj  +  dxSi  +  -jir  duxSi 


(31) 


The  source  term  discretization  is,  therefore,  not  conservative,  which 
is  natural  and  appropriate  because  the  source  term  does  not  have  a 
conservative  property  and  should  not  be  discretized  in  such  a  way. 
If  it  is  conservative,  the  global  sum  of  the  cell-residual  for  the  source 
term  will  depend  only  on  the  boundary  data  (telescoping  property), 
which  is  wrong  for  the  source  term. 

We  now  show  that  the  truncation  error  (T.E.)  of  the  resulting  RD 
scheme  with  the  above  divergence  formulation  of  the  source  term 
(with  m  =  3)  is  in  fact  third-order.  Again,  for  demonstration 
purposes  we  consider  the  second  equation  (i.e.  dzp);  the  same  pro¬ 
cess  can  be  repeated  for  the  first  equation.  We  first  expand  the 
source  fluxes  around  the  node  i: 

jf "  dxfsdx=fl,  -jf  =  hRdJs  +  +  0(h4), 

(32) 


- o - o - o - 

i-1  hL  i  hR  i+1 

Fig.  3.  Cell  residual  evaluation  with  divergence  form  of  the  source  term. 


where  the  d*fs  is  given  by  Eq.  (28),  and 

djs  =  dxS  +  h-j-  dmS  +  ^  d(x4)S  +  0(h4) ,  (33) 

Ij3 

<WS  =  <9x*S  +  hRdmS  +  h2Rdl4)S  +  g  df  >s  +  0(h4) .  (34) 


Using  Eq.  (29),  we  evaluate  the  truncation  error  of  the  equation  dzp 
by  substituting  Eqs.  (33)  and  (34)  into  Eq.  (32)  and  expanding  all  the 
terms  in  Eq.  (29)  around  the  node  i: 


T.E.(dTp)  = 


.i2 

(dxlli  —  Pi)  +  2  ( —  dxPi)  +  g  [dxxxUi  —  <)xxPi) 


Fr 


+ 


^(d^ui-\4dxxxPi)  +  0(h4) 


/Tr 


=  0  +  0(h3), 


(35) 


where  the  first  three  terms  vanish  as  they  satisfy  our  original 
equation.  The  presence  of  the  last  term  confirms  that  the  proposed 
divergence  formula  with  m  =  3  makes  our  scheme  third-order 
accurate.  The  same  result  is  obtained  for  the  first  equation  with 
the  divergence  formulation  of  the  corresponding  source  terms, 
and  therefore  the  process  is  not  repeated  here. 

We  remark  that  the  cost  of  this  new  third-order  accurate  RD 
hyperbolic  advection-diffusion  scheme  using  the  divergence 
formulation  of  the  source  terms  is  only  the  evaluation  of  the  first 
and  second  derivatives  of  the  source  terms;  i.e.,  dxS  and  c^S.  Note 
that  the  source  flux  is  third-order  accurate  (see  Eq.  (28)).  Thus, 
according  to  Eqs.  (30)  and  (31),  we  need  second-order  and  first- 
order  accurate  discretization  for  dxS  and  respectively.  We 
derive  these  discretization  by  fitting  a  quadratic  function  through 
the  node  i  and  its  two  neighbors  i  -  1  and  i  +  1  to  arrive  (after  some 
algebra)  at  the  following  formulas  that  are  applicable  to  the  inter¬ 
nal  nodes  of  general  arbitrary  grids  (uniform  and  nonuniform): 


dxSi  =  +  hl)Si  +  h2LSi+ 1  +  0(h2^ 


dxxSi  = 


hRhL(hR  +  hL) 

hRSi-i  -  (hR  +  hjjSi  +  hiSj+i 

hRhL(hR  +  hL)/2 


+  0(h), 


(36) 

(37) 


where  hR  and  hL  are  defined  in  Eq.  (10).  The  corresponding  one¬ 
sided  formulas  for  the  boundary  nodes  of  general  arbitrary  grids  are 


r,  c  — hr* (far  +  2hR)Sj  +  (hR  +  hL-)2S2  —  hpS3 

C'xm  = - u  u  7u —  +  ' - r  ), 


dxSN  = 
<9*xSi  = 
&xxSn  — 


hRhL*{hR  +  hL.) 
hi.{hi  +  2hR*)SN  —  ( hR •  +  hL)2SN_j  +  hp.SN_2 
hR'hL(hR>  +  hL) 
hL'Sj  -  (hR  +  hL>  )S2  +  hRS3 
hRhL<(hR  +  hr)/2 
—h/rSn  +  (/tr*  +  hj.)Siv_i  —  h/Siv- 2 


hR<hL(hR’  +  hL)/2 


(38) 

+  0(h2),  (39) 

+  0(h),  (40) 

0(h),  (41) 


where  hr  and  hR.  are  defined  as 

hL>  =  x3  -  x2 ,  hR.  =  xN_,  -  xN_2 ■  (42) 

It  is  clear  from  the  above  formulas  that  each  derivative  can  be 
computed  in  a  three-point  stencil.  Consequently,  the  residual  at  a 
node  is  defined  in  a  five-point  stencil  in  the  interior,  a  four-point 
stencil  at  the  nodes  adjacent  to  the  boundary,  and  a  three-point 
stencil  at  the  boundary  nodes.  In  the  next  section,  we  demonstrate 
that  a  fourth-order  scheme  can  be  constructed  without  extending 
the  stencil. 
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4.2.  Fourth-order  scheme  with  divergence  form 


Here,  we  show  how  a  simple  modification  to  the  presented 
divergence  form  upgrades  the  scheme  order  by  one  order;  that  is 
our  third-order  scheme  becomes  fourth-order  with  no  additional 
cost.  To  gain  an  order,  we  propose  the  divergence  form  of  the 
source  terms  with 

m> 3  /_-]  \ n— 1 

fs  =  Y,(-^r-  (x-x)n^-'S,  (43) 

n= 1 

=  (x  -  x)S  -  i  (x  -  x)2dxS  +  g  (x  -  (x  -  x)4dxxxS  +  ..., 

(44) 

where  x  —  (x,-  +  xi+i) /2.  Note  that  the  previous  divergence  form  was 
formed  around  the  x ,  while  this  new  formulation  defines  the  diver¬ 
gence  function  around  the  mid-point  x.  Similar  to  the  divergence 
form  presented  in  the  previous  section,  the  source  term  and  its 
derivatives  are  not  constants  but  functions  that  vary  in  space.  And 
therefore,  when  the  source  flux  is  differentiated  it  recovers  the  ori¬ 
ginal  source  term  up  to  0(hm )  around  the  node  i.  For  m  =  3,  we  have: 

dxfs  =  S  +  ^^-dxxxS  =  S  +  0(h3).  (45) 

Again,  it  does  not  limit  the  maximum  order  of  accuracy  for  numer¬ 
ical  schemes,  and  it  is  possible  to  construct  even  higher-order 
schemes. 

We  now  show  that  this  third-order  source  divergence  formula¬ 
tion  will  result  in  a  fourth-order  accurate  RD  scheme.  We  do  this 
similarly  to  the  process  explained  in  the  previous  section  except 
that  the  second  and  third  derivatives  of  the  fs  are  now  defined  as: 

djs  =  dxS  +  ^y^-  dmS  +  9<4»S  +  0(h4) ,  (46) 


d,oJs  =  d^S  +{x-  x)dmS  +  (x-  x)2d{4)S  +  (X  gX)  45)S  +  0(h4). 


(47) 


Again,  for  discussion  purposes,  we  consider  the  second  equation  of 
the  hyperbolic  advection-diffusion  system,  dzp.  The  truncation 
error  of  the  p  equation  after  substituting  Eqs.  (46)  and  (47)  into 
Eq.  (32)  becomes: 


T.E.(dtp)  = 


( dxUi  - +  y  {dxxUi  -  dxPi)  +  y  (dxxxUi  -  dxxPi ) 


+ 


^(d^Ui  -  d^pj  +  ^dfui  -  ^<4>Pi)  +  0(h5 


=  0  +  0(fr), 


/T  r 

(48) 


which  is  fourth-order  accurate  because  of  the  cancellation  of  the 
first  four  terms  of  the  truncation  error  equation.  Another  great 
property  of  this  new  divergence  formulation  is  the  equal  evaluation 
of  the  /))  |  -ff  regardless  of  whether  the  source  flux  is  being 
transferred  to  the  node  i  or  i+1.  This  greatly  simplifies  the 
implementation  of  this  form  of  divergence  formulation. 

Note  that  the  identical  residual  is  distributed  to  the  left  and 
right  nodes  in  this  scheme.  It  then  appears  that  the  source  term 
discretization  is  conservative.  However,  it  is  actually  not  conserva¬ 
tive  for  the  source  term  because  the  cell-residual  for  the  source 
term  depends  on  the  coordinate  of  the  mid-point  of  the  cell  and 
thus  it  is  not  telescoping.  Again,  this  is  natural  and  appropriate 
for  the  source  term,  which  has  no  conservation  property  and 
should  be  local. 


5.  Higher-order  RD  scheme  with  generalized  trapezoidal  rule 
(third-,  fourth-  and  sixth-order) 


In  the  previous  section,  we  reformulated  the  original  hyperbolic 
advection-diffusion  system  with  a  generalized  divergence  form  of 
the  source  terms  and  arrived  at  third-  and  fourth-order  RD 
schemes.  We  specifically  showed  that  both  the  third-  and  the 
fourth-order  RD  schemes  developed  with  the  divergence 
formulation  of  the  source  terms  require  the  evaluation  of  the  first 
and  second  derivatives  of  the  source  terms.  Fifth-  and  higher-order 
schemes  can  be  systematically  constructed  with  the  divergence  for¬ 
mulation,  but  may  require  the  evaluation  of  third-  and  higher-order 
derivatives  that  would  extend  the  stencil  to  the  neighbors  of  the 
neighbors  and  beyond.  From  a  practical  point  of  view,  such  high- 
order  schemes  are  not  very  attractive,  and  thus  not  considered  here. 

In  this  section,  we  introduce  a  different  technique  to  develop 
even  higher-order  schemes  without  the  need  to  evaluate  gradients 
beyond  the  second-derivatives.  We  also  show  that  with  this  new 
technique  the  fourth-order  RD  scheme  is  even  less  computation¬ 
ally  intensive  than  the  third-  and  fourth-order  RD  schemes  that 
are  developed  with  the  divergence  formulation  of  the  source 
terms. 

Consider  the  vector  form  of  our  first-order  hyperbolic  advec¬ 
tion-diffusion  system: 


_l  /Y—  =  s 

dr +  dx 


(49) 


We  showed  in  Section  3  that  the  source  term  discretization  with  the 
trapezoidal  rule  provides  a  second-order  accurate  scheme.  This 
trapezoidal  rule  can  be  written  as 

J‘+'  Sdx  =  y  (Sj,  +  Sr),  (50) 


where  for  the  second-order  scheme  SL  =  S,  and  SR  =  Si+1 ;  i.e.,  the 
arithmetic  averaging  of  the  source  terms  between  the  left  and  the 
right  nodes  (Fig.  4).  Generalizing  the  trapezoidal  rule,  we  propose 
the  following  formula  for  the  left  and  right  source  terms,  SL  and 
SR  respectively: 

SL  =  St  +  C\dxSt  +  C^Sj,  (51) 

SR  =  SM  +  C*dxSi+ 1  +  C%Si+1 ,  (52) 

where  and  are  constants  of  0(h),  and  CL2  and  C2  are  of  0(h2).  A 
somewhat  similar  approach  is  taken  in  Ref.  [15],  introduced  for 
upgrading  the  order  of  accuracy  of  a  finite-volume  scheme,  where 
not  only  the  source  term  but  also  the  interface  flux  need  to  be  mod¬ 
ified  for  high-order  accuracy.  We  find  these  constants  by  making 
sure  that  the  nodal  residuals  of  the  two  equations  in  our  hyperbolic 
system  are  accurate  up  to  the  desired  order  of  accuracy  of  the 
scheme  we  are  seeking  to  develop.  With  the  above  equation,  the 
maximum  possible  order  of  accuracy  is  six.  We  also  note  that  the 
fifth-order  RD  scheme  with  the  above  proposed  source  term  dis¬ 
cretization  becomes  mathematically  sixth-order. 


5.1.  Third-order  RD  scheme  with  generalized  trapezoidal  rule 


A  third-order  RD  scheme  can  be  obtained  if  the  coefficients  of 
the  proposed  discretization  formulation  satisfy  the  following 
equations: 


-O- 


i-1 


Fig.  4.  Source  term  discretization. 


-o 


i-1 
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c;  +  ct  =  0,  CX  +  C2  +  C*  =  -^,  CfflR  +  2Cf^-ft,  (53) 

o  o 


which  can  be  derived  from  the  truncation  error  analysis.  Here,  there 
are  four  unknowns  for  three  equations;  thus  there  are  infinite  com¬ 
binations  with  the  Cj  as  a  free  parameter.  We  also  note  that  there  is 
no  relation  between  this  third-order  RD  scheme  and  the  third-order 
scheme  developed  with  the  divergence  formulation;  the  nodal 
residual  of  the  source  terms  obtained  with  the  previous  third-order 
scheme  depends  only  on  the  information  from  its  immediate  node 
neighbor  and  not  the  node  itself.  On  the  other  hand,  the  new 
third-order  RD  scheme  that  is  developed  with  a  correction  to  the 
trapezoidal  rule  depends  both  on  its  immediate  node  neighbor 
and  the  node  itself.  Consequently,  it  is  not  possible  to  reproduce 
the  previous  third-order  scheme  by  any  choice  of  the  coefficients. 
Note,  however,  that  this  scheme  has  the  same  stencil  with  the 
previous  third-order  scheme:  a  five-point  stencil  in  the  interior,  a 
four-point  stencil  at  the  nodes  adjacent  to  the  boundary,  and  a 
three-point  stencil  at  the  boundary  nodes.  The  same  is  true  for 
the  fourth-order  scheme  derived  in  the  next  section. 

Assuming  =  -hR/ 6,  the  left  and  right  source  functions  of  the 
source  discretization  may  be  chosen  as: 


SR  =  S!+1 


(54) 

(55) 


We  can  now  evaluate  the  truncation  error  of  the  hyperbolic  advec- 
tion-diffusion  system.  Here  we  show  the  truncation  error  of  the 
second  equation;  the  same  will  be  true  for  the  first  equation  as  well. 
The  cell  residual  of  the  second  equation  with  the  third-order  RD 
scheme  with  the  generalized  trapezoidal  rule  is: 

<  =  -a(uM  -  ut)  +  v(pM  -  p,-)  +  y  (SL  +  SR), 

hR 

=  -a(ui+1  -  u,)  +  v(Pi+1  -  p,)  +  y  (S,-+1  +  Si) 

-  ft  (&si+1  -  dxSi)  -  ^ -  (aMSi+1  -  dxxSi).  (56) 


We  then  expand  the  cell  residual  around  the  node  i  to  arrive  at  the 
following  truncation  error: 

T.E.(dxUi)  =  ( -adM  +  vdxpt  +  S,)  +  y  (-ac^U;  +  Vc^p,-  +  dxSt) 

+  q  (—  QdxxxUi  +  VdioaPj  +  dxxSi) 

h3 

+  2^  (-acfti,  +  Vd^Pi  +  0.988  dmSt) 

+  0(h4)  =  0  +  0(h3),  (57) 


where  the  first  three  terms  vanish  but  the  last  term  makes  the 
scheme  third-order  accurate.  As  we  will  show  next,  smaller  values 
for  the  coefficients  C2  or  C2  move  the  schemes  to  fourth-order  accu¬ 
rate.  This  is  further  shown  in  the  following  section. 


Ci=+ft,  C‘=0,  C?  =  -ft,  C*  =  0.  (59) 

As  will  be  shown  shortly,  it  leads  to  a  fourth-order  RD  scheme  with 
the  source  term  discretization  by  Eqs.  (51)  and  (52): 

SL  =  Si  +  h*dxSt,  (60) 

SR  =  SM-^dxSi+u  (61) 

It  is  remarkable  that  the  fourth-order  scheme  does  not  require  the 
evaluation  of  the  second  derivatives  and  is  thus  less  expensive  than 
the  third-order  schemes  developed  in  the  previous  sections.  It  is,  of 
course,  possible  to  develop  a  fourth-order  RD  scheme  with  addition 
of  the  second  derivatives  in  the  source  term  discretization.  For 
example  the  following  constants  will  also  result  in  a  fourth-order 
RD  scheme: 


rr  _  hR  L  _  Iir  r  _  hR  R  _  hR 
C'“+T’  L2~+24  4~’  L2_+24' 


(62) 


which  makes  this  fourth-order  RD  scheme  identical  to  the  fourth- 
order  RD  scheme  constructed  by  the  divergence  formulation  in  Sec¬ 
tion  4.2.  Although  it  is  possible  to  construct,  these  fourth-order 
schemes  are  not  very  attractive  as  they  require  second  derivatives 
of  the  source  term,  and  therefore  not  considered  in  this  paper. 

We  now  prove  the  fourth-order  accuracy  of  the  scheme  by  eval¬ 
uating  the  cell  residual  of,  for  example,  the  first  equation;  i.e. 

<  =  -a(uM  -  u,)  +  v(pM  -  Pi)  +  y  (Si  +  Sr), 

=  -a(uM  -  Ui)  +  v(pi+1  -  Pi)  +  y  (Si+1  +  Si)  -  ft  (0*Si+1  -  dxSi). 

(63) 


Expanding  the  cell  residual  around  the  node  i,  we  obtain  the  trun¬ 
cation  error  of  the  first  equation  (after  some  algebra)  as 

T  ,E.(drUi)  =  (—adxUj  +  vdxPi  +  S,)  +  —  (—cii)xxUi  +  vdxxpl  +  dxS\) 

+  y  {~ddxxxUi  +  VdioaPj  +  0XxS,) 

+  2^(~ad^Ui  +  vd^Pj  +  dxxxSj ) 

+  U0<'~a9^Ui  +  Vd*)pi  +  g^>S') 

+  0(h5)  =  0  +  0(h4).  (64) 


where  the  first  four  terms  of  the  above  equation  vanish  (because  of 
consistency  relations).  Similar  conclusion  is  obtained  for  the  second 
equation  and  therefore  it  is  not  repeated  here.  Also,  note  that  the 
accuracy  is  achieved  for  general  arbitrary  grids.  The  additional  cost 
of  upgrading  the  second-order  scheme  to  the  fourth-order  RD 
scheme  is  only  due  to  the  evaluation  of  the  second-order  accurate 
first  derivative  of  the  source  terms.  The  general  second-order  accu¬ 
rate  first  derivative  formulation  for  arbitrary  grids  is  provided  in 
Section  3.1. 


5.2.  Fourth-order  RD  scheme  with  generalized  trapezoidal  rule 


Fourth-order  accuracy  is  achievable  by  a  simple  adjustment  to 
the  condition  given  in  Eq.  (53): 


ci  +  c?  =  o,  cX  +  c‘  +  c«  =  -|, 


CR 
L1 1 


rR  _  llR 

—  T2 


(58) 


Again,  the  coefficients  cannot  be  determined  uniquely,  and  there 
are  many  solutions.  A  particularly  interesting  solution  is  the 
following: 


5.3.  Sixth-order  RD  scheme  with  generalized  trapezoidal  rule 

In  this  section,  we  present  a  new  sixth-order  RD  scheme  with 
relatively  similar  cost  to  our  newly  introduced  third-order  RD 
schemes.  Specifically,  we  show  that  there  exists  a  unique  fifth- 
order  RD  scheme  with  the  generalized  trapezoidal  rule,  and  that 
mathematically  becomes  sixth-order  RD  scheme.  Seeking  fifth- 
order  accuracy  from  the  fourth-order  RD  scheme  given  in  Eq. 
(58),  we  find  an  additional  constraint: 


138 


A  Mazaheri,  H.  Nishikawa / Computers  &  Fluids  102  (2014)  131-147 


r'R  _  '  lR 

C2~~20' 


(65) 


Then,  the  number  of  constraints  (four)  matches  the  number  of 
unknown  coefficients  (four).  In  this  case,  therefore,  there  exists  a 
unique  solution: 


(66) 


Interestingly,  the  above  unique  coefficients  satisfy  the  following 
constraint  as  well,  which  is  the  requirement  for  sixth-order 
accuracy: 


i  h  i  rR  —  R 

4  hR  +  c2--30- 


(67) 


Thus,  we  have  developed  a  sixth-order  RD  scheme  with  the  general¬ 
ized  trapezoidal  rule.  We  remark  that  these  coefficients,  unlike  the 
ones  proposed  for  the  third-order  and  fourth-order  schemes,  are 
unique.  With  the  proposed  constants  for  the  sixth-order  RD  scheme, 
the  source  term  evaluations  at  the  nodes  i  and  i  +  1  become: 


SL  =  Si+h^dxSi  +  ^dxxSi,  (68) 

SR  =  SM  -  h dxSM  +  ^ aMSi+1 .  (69) 

Note  that  it  requires  only  the  first  and  second  derivatives  of  the 
source  term,  and  no  higher-order  derivatives  are  required. 

Similar  to  the  process  explained  for  the  fourth-order  scheme, 
we  prove  the  order  of  accuracy  of  the  scheme  using  the  above  pro¬ 
posed  source  discretization  by  first  evaluating  the  cell  residual  of 
the  hyperbolic  system.  Consider  the  cell  residual  of,  for  example, 

the  first  equation  (@‘u)  • 

<  =  -a(ui+1  -  Ui)  +  V(Pi.l  ~  Pi)  +  y  (St  +  SR), 

fjn 

=  -a(uM  -  ui)  +  v(pi+1  -  p^  +  y  (Si+1  +  St) 

“  ^  (dxsM  -  dxSi)  +  ^  (3„sit.,  +  dxxSi) .  (70) 

Expanding  the  cell  residual  around  the  node  i,  we  obtain  the  trun¬ 
cation  error  of  the  first  equation  (after  some  algebra)  as 


T.E.(dzUi)  =  {-adM  +  vdxPi  +  S,) 

fin 

+  ~2  ( —ud^Ui  +  vdxxPi  +  dxSi) 

+  q  (  -  UiJxxxU,  +  VdxxxPj  +  dmSj) 

+  ^  (-acf  Ui  +  vd^  +  dmS,) 

+  4  (-ad?Ui  +  vdVPi  +  d?Si) 

+  *L(-at$>Ui  +  V#?pi  +  &St) 

+  5M0  (‘a^)Ui  +  Vd*)pi  +  S^S')  +  °(h7) 

=  0  +  0(h6).  (71) 


A  similar  result  is  obtained  for  the  second  equation,  and  is  thus 
omitted  here,  for  brevity.  The  above  truncation  error  analysis 
reveals  that  this  powerful  sixth-order  RD  scheme  could,  in  practice, 
produce  almost  seventh-order  accurate  results,  if  the  sixth  deriva¬ 
tive  of  the  source  term  becomes  small;  the  sixth-order  term  is  only 
due  to  the  presence  of  (h^/100800)9*6)Sj  in  the  last  term  of  Eq.  (71). 
We  will  show  such  interesting  results  in  the  following  section. 


We  emphasize  that  the  sixth-order  RD  scheme,  similar  to  the 
third-order  RD  scheme,  requires  the  evaluation  of  the  first  and  sec¬ 
ond  derivatives  of  the  source  terms.  However,  these  derivatives  are 
now  required  to  be  evaluated  with  higher-order  accuracy.  For  this 
sixth-order  RD  scheme,  we  need  third-order  and  second-order 
accurate  evaluations  of  the  first  and  the  second  derivatives  of  the 
source  terms,  respectively,  which  can  be  obtained  by  a  cubic  fit. 
This  makes  the  developed  sixth-order  scheme  slightly  more  expen¬ 
sive  than  the  third-order  scheme.  Nevertheless,  the  proposed 
sixth-order  scheme  is  exceptionally  simple  and  affordable. 

6.  Results 

In  this  section  we  present  the  results  in  three  categories:  (1) 
steady  advection-diffusion  equation  for  high  Reynolds  (or  Peclet) 
number  applications,  (2)  unsteady  linear  advection-diffusion, 
and  (3)  unsteady  nonlinear  advection-diffusion  problems.  We 
obtain  the  results  with  all  the  proposed  RD  schemes  and  compare 
them  with  the  second-order  RD  scheme.  The  order  of  accuracy 
results  are  also  compared  and  presented  for  each  example. 

For  all  cases,  we  employ  the  Newton-GS  solver  as  described  in 
Section  3.  It  is  essentially  Newton’s  method  for  the  second-order 
scheme,  but  an  approximate  Newton  method  for  higher-order 
schemes  because  the  Jacobian  matrix  derived  from  the  second- 
order  scheme  is  used  for  all  higher-order  schemes.  For  unsteady 
problems,  the  same  solver  is  used  to  solve  the  implicit-residual 
equations  or  equivalently  to  compute  the  pseudo  steady  solution 
at  every  physical  time  step. 

6.1.  Steady  linear  advection-diffusion 

Consider  the  advection-diffusion  equation  in  x  e  (0, 1)  with 
u(0)  =  0  and  u(l)  =  0: 

dtu  +  adxu=vdxxu+s(x),  (72) 

where 

s(x)  =  vn2  sin(7ix)  +  a7tcos(7ix).  (73) 

The  above  problem  has  a  fixed  analytical  solution  of  sin(7tx)  for  any 
advection  and  diffusion  coefficients.  We  solved  this  problem  using 
the  hyperbolic  advection-diffusion  method  discussed  in  Section  2 
with  the  proposed  high-order  RD  schemes  (see  Fig.  5.  Note  that 
we  intentionally  used  a  very  coarse  grid  to  show  that  the  high-order 
schemes  could  produce  a  much  more  accurate  solution  on  such  a 
“bad"  grid.)  We  used  ranges  of  nonuniform  grids  and  solved  the 
problem  with  the  Newton-GS  method.  For  each  Newton  iteration, 
the  GS  relaxation  were  conducted  until  three  orders  of  magnitudes 
reduction  is  achieved  for  the  linear  system.  The  computations  were 
continued  until  the  residuals  of  both  the  solution  u  and  the  solution 
gradient  p  were  reduced  by  eight  orders  of  magnitude.  The  conver¬ 
gence  results  are  shown  in  Table  2,  where  the  convergence  of  the  GS 
relaxation  is  clearly  0(N)  and  is  independent  of  the  order  of  accu¬ 
racy  of  the  RD  scheme.  Also,  the  solution  was  converged  with  the 
same  number  of  GS  relaxations  and  Newton  iterations  regardless 
of  the  RD  scheme  order.  Note  that  the  solutions  were  converged 
with  a  very  small  number  of  Newton  iterations,  typically  less  than 
ten;  this  is  exceptionally  remarkable  for  the  approximate  Jacobian 
(second-order)  formulation  for  higher-order  schemes. 

The  accuracy  of  the  proposed  RD  schemes  were  verified  by  com¬ 
puting  the  Li  =  (Ur1  ~  Ui)  /N-  These  results  are  shown  in 
Fig.  6  for  the  third-,  fourth-,  and  sixth-order  hyperbolic  RD 
schemes  (see  also  Tables  3-5).  Note  that  when  the  differences 
between  the  numerical  results  and  exact  values  approach  the 
machine  zero,  the  LI  slope  approaches  zeroth-order;  these  results 
are  thus  omitted  from  the  Table  5.  The  RD  schemes  developed  with 
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(a)  u{x)  (b)  p(x )  =  du/dx 


Fig.  5.  Comparisons  between  the  exact  and  the  numerical  results  using  the  proposed  high-order  RD  schemes  for  the  steady  linear  advection-diffusion  problem  (Re  =  1)  on  a 
nonuniform  grid  with  N  =  10. 


Table  2 

Steady  linear  advection-diffusion  problem,  Re  =  1  (Convergence  criteria:  residuals  <  1CT8.) 


Number  of  nodes 

RD  scheme  order 

GS  relaxations/Newton  iteration 

High-order  technique 

Newton  iteration 
High-order  technique 

RD-D 

RD-GT 

RD-D 

RD-GT 

50 

3rd 

168 

169 

10 

10 

4th 

168 

168 

10 

10 

6th 

- 

168 

- 

10 

100 

3rd 

325 

325 

8 

8 

4th 

325 

325 

8 

8 

6th 

- 

325 

- 

8 

200 

3rd 

670 

670 

7 

7 

4th 

670 

670 

7 

7 

6th 

- 

670 

- 

7 

300 

3rd 

1015 

1015 

7 

7 

4th 

1015 

1015 

7 

7 

6th 

- 

1015 

- 

7 

500 

3rd 

1703 

1703 

7 

7 

4th 

1703 

1703 

7 

7 

6th 

- 

1703 

- 

7 

1000 

6rd 

3416 

3416 

7 

7 

4th 

3416 

3416 

7 

7 

6th 

- 

3416 

- 

7 

Third-Order  RD  Schemes  Fourth-Order  RD  Schemes  Sixth-Order  RD  Scheme 


log10h  log10h  logI0h 

(a)  third-order  (b)  fourth-order  (c)  sixth-order 


Fig.  6.  L,  error  of  the  proposed  high-order  RD  hyperbolic  schemes  for  the  steady  linear  advection-diffusion  problem. 
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the  divergence  formulation  are  denoted  by  RD-D,  while  the 
schemes  developed  with  the  generalized  trapezoidal  rule  are 
denoted  with  RD-GT.  The  results  show  that  all  the  RD  schemes 
achieve  the  design  order  of  accuracy  for  both  the  solution  u  and 
the  solution  gradient  p.  It  can  be  seen  also  that  the  difference 
between  the  two  versions  of  third-order  RD  schemes  is  minor, 
and  the  same  is  true  for  the  fourth-order  RD  schemes. 

6.2.  Steady  boundary  layer  problem 

Consider  the  advection-diffusion  equation  in  x  e  (0, 1)  with 
u( 0)  =  0  and  u(l)  =  1 : 

d  tu  +  adxu  =  vdxxii  +  s(x),  (74) 

where 

71 

s(x)  =  —  [ac°s(7ix)  +  7tvsin(7tx)],  Re  =  a/v.  (75) 

This  is  a  boundary  layer  problem  with  a  non-trivial  steady  state 
solution  in  the  diffusion  limit  as  a  result  of  the  source  term  addition 
[2],  This  equation  develops  a  very  narrow  boundary  layer  near  the 
right  boundaiy  (x  =  1 )  when  the  advection  term  becomes  domi¬ 
nant.  The  exact  steady  state  solution  to  this  problem  is  given  by 
(see  also  Fig.  7). 

p-Re  _  p(x-l)Re  i 

“““(*)  =  e_fe_1  +^sin(7tx).  (76) 

We  chose  various  Re  values  ranging  from  1  to  106  and  solved 
the  equation  on  nonuniform  grid  sizes  up  to  100,000  nodes.  Like 
the  previous  problem,  the  solutions  were  obtained  with  the 
Newton-GS  method.  Within  each  Newton  iteration,  the  GS  relaxa¬ 
tion  were  conducted  until  three  orders  of  magnitude  reduction  is 
achieved  for  the  linear  system,  and  the  residuals  of  both  the  solu¬ 
tion  and  the  solution  gradient  were  reduced  by  eight  orders  of 
magnitude.  The  convergence  data  are  given  in  Table  6  for  the  pro¬ 
posed  high-order  RD  schemes.  The  data  are  shown  for  the  high- 
order  schemes  developed  with  the  generalized  trapezoidal  rule. 
Similar  results  were  obtained  with  the  RD  schemes  developed  with 
the  divergence  formulation,  and  therefore  not  shown.  But  in  gen¬ 
eral,  the  high-order  RD  schemes  with  the  generalized  trapezoidal 
rule  perform  better  than  the  ones  developed  with  the  divergence 
formulation.  The  latter  third-order  RD-D  scheme  not  only  produce 
slightly  larger  errors  as  shown  in  Fig.  6,  but  also  encounter  some 
convergence  difficulties  (particularly  with  time-dependent 
problems);  the  convergence  problem  seems  originated  from  the 
lack  of  diagonal  dominance  caused  by  the  particular  structure  of 
the  source  term  discretization  as  shown  in  Section  4.1.  We  will 
therefore  mainly  focus  on  the  RD-GT  schemes  for  unsteady  cases. 
The  O(N)  linear  dependency  of  the  GS  relaxations  on  the  grid  size 
is  also  demonstrated,  which  is  a  consequence  of  solving  the  advec¬ 
tion-diffusion  equation  as  a  hyperbolic  system.  We  emphasize  that 
this  is  remarkable  because  the  linear  convergence  is  retained  for 
any  irregular  grid  in  any  dimensions  (N  is  approximately  the 
number  of  nodes  in  each  coordinate  direction  in  two  and  three 


Table  3 

Spatial  accuracy  with  the  third-order  RD-GT  scheme  for  the  steady  linear  advection- 
diffusion  problem  (Re  =1). 


Number  of  nodes 

Li  error  of  u 

Order 

Li  error  of  p 

Order 

10 

9.54E-03 

8.59E-02 

25 

4.63E-04 

3.30 

4.41  E-03 

3.24 

50 

4.85E-05 

3.25 

4.63E-04 

3.25 

100 

5.46E-06 

3.15 

5.16E-05 

3.17 

200 

6.46E-07 

3.08 

6.04E-06 

3.09 

300 

1.88E-07 

3.04 

1.75E-06 

3.05 

500 

3.99E-08 

3.03 

3.71E-07 

3.04 

1000 

4.93E-09 

3.02 

4.57E-08 

3.02 

Table  4 

Spatial  accuracy  with  the  fourth-order  RD-D  scheme  for  the  steady  linear  advection- 
diffusion  problem.  Re  =  1.  (Note:  fourth-order  RD-D  and  RD-GT  schemes  produce 
almost  identical  result.) 


Number  of  nodes 

Li  error  of  u 

Order 

Li  error  of  p 

Order 

10 

8.38E-03 

6.36E-02 

25 

2.84E-04 

3.69 

1.82E-03 

3.88 

50 

1.74E-05 

4.03 

1.02E-04 

4.16 

100 

1.04E-06 

4.07 

5.77E-06 

4.15 

200 

6.23E-08 

4.06 

3.35E-07 

4.11 

300 

1.21E-08 

4.04 

6.44E-08 

4.07 

500 

1.55E-09 

4.02 

8.15E-09 

4.05 

1000 

9.30E-1 1 

4.06 

4.94E-10 

4.04 

Table  5 

Spatial  accuracy  with  the  sixth-order  RD-GT  scheme  for  the  steady  linear  advection- 
diffusion  problem  (Re  =  1 ). 


Number  of  nodes 

It  error  of  u 

Order 

It  error  of  p 

Order 

10 

2.60E-03 

1.75E-02 

25 

3.34E-05 

4.75 

1.98E-04 

4.89 

50 

4.87E-07 

6.10 

2.87E-06 

6.11 

100 

4.32E-09 

6.82 

2.65E-08 

6.76 

200 

3.09E-11 

7.13 

1.78E-10 

7.22 

300 

5.80E-12 

- 

3.39E-11 

- 

500 

3.60E-12 

- 

1.97E-11 

- 

dimensions)  as  demonstrated  in  Refs.  [1-3, 6, 4, 5].  It  leads  to  orders 
of  magnitude  faster  convergence  in  comparison  with  conventional 
methods  whose  convergence  is  typically  0(N2)  as  discussed  also  in 
the  previous  paper  [6],  Furthermore,  the  proposed  high-order  RD 
schemes  are  extremely  efficient  as  the  solutions  were  obtained 
with  only  a  small  number  of  Newton  iterations:  less  than  10 
iterations  to  reduce  the  residual  by  eight  orders  of  magnitude. 
Moreover,  the  number  of  GS  relaxations  and  Newton  iterations 
are  essentially  independent  of  the  scheme  order.  Considering  the 
fact  that  the  cost  of  one  GS  relaxation  is  significantly  cheaper  than 
one  Newton  iteration,  we  find  that  the  developed  high-order  RD 
schemes  are  extremely  powerful  and  efficient.  Finally,  as  in  the 
previous  work  [6],  we  remark  that  the  high-Re  cases  required 
extremely  fine  grids  to  meet  the  well-known  requirement  on  the 
mesh  Reynolds-number  [2].  If  desired,  the  computations  can  be 
performed  on  substantially  coarser  grids  with  more  aggressive 
and  customized  grid  stretching.  However,  we  simply  refined  the 
grid  to  meet  the  mesh  Reynolds-number  requirement  because 
our  method  is  powerful  enough  to  solve  the  problem  very 
efficiently  (i.e.,  less  than  10  Newton  iterations)  even  on  such  dense 
grids  for  all  third-,  fourth-,  and  sixth-order  schemes.  The  ability  to 
efficiently  solve  the  problem  on  highly  refined  grids  is  a  great 
advantage  of  these  schemes. 

The  order  of  accuracy  of  the  proposed  RD  schemes  were  also 
verified  for  this  problem.  Figs.  8  shows  the  L i  error  convergence 
results,  where  h  is  the  representative  mesh  spacing  defined  by 
h  =  /(N-  1).  For  discussion  purposes,  we  present  the  accuracy 
plots  for  Re  =  1  and  Re  =  10;  similar  results  were  obtained  for 
other  Reynolds  numbers.  These  results  verify  the  order  of  accuracy 
of  the  proposed  high-order  RD  schemes  (i.e.,  third-,  fourth-,  and 
sixth-order)  for  all  the  variables  and  the  gradients  at  all  the  grid 
nodes  including  the  boundary  nodes  (see  also  Tables  7-9).  Note 
that  when  the  differences  between  the  numerical  results  and  exact 
values  approach  the  machine  zero,  the  L i  slope  approaches 
zeroth-order;  these  data  are  omitted  from  the  Tables  8  and  9. 

6.3.  Unsteady  linear  advection-diffusion 

Consider  the  time-dependent  advection-diffusion  equation  in 
xe  (0.1) 
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Re=10 


0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1 

X 


(b)  p(x)  =  du/dx 


Fig.  7.  Comparisons  between  the  exact  and  the  numerical  results  using  the  proposed  high-order  RD  schemes  for  the  steady  boundary  layer  problem  (Re  =  10)  on  a 
nonuniform  grid  with  N  -  10. 


Table  6 

Steady  boundary  layer  problem  (Convergence  criteria:  Residuals  <  1(TS.) 


log10Re 

Number  of 
nodes 

RD-GT 

Scheme  Order 

GS  relaxations/ 
Newton  iteration 

Newton 

iteration 

0 

50 

3rd 

163 

8 

4th 

163 

8 

6th 

0 

100 

3rd 

324 

7 

4th 

324 

7 

6th 

0 

500 

3rd 

1647 

7 

4th 

1647 

7 

6th 

1 

100 

3rd 

178 

7 

4th 

178 

7 

6th 

2 

100 

3rd 

44 

7 

4th 

44 

7 

6th 

3 

500 

3rd 

42 

8 

4th 

43 

7 

6th 

4 

1000 

3rd 

18 

9 

4th 

19 

7 

6th 

5 

10,000 

3rd 

24 

9 

4th 

21 

7 

6th 

6 

100,000 

3rd 

18 

9 

4th 

19 

7 

6th 

dtu  +  adxii  =  vdxxii.  (77) 

The  above  equation  with  the  following  initial  condition: 

u(x,  t  =  0)  =  sin(icx),  (78) 

where  k  is  an  arbitrary  constant,  has  the  following  exact  solution 
with  a  periodic  boundary  condition: 

uexact(x,  t)  =  e^vt  sin(K(x  -  at)).  (79) 

A  non-periodic  solution  also  exists  for  the  following  oscillatory 
boundary  conditions 

u(0,  t)  =0,  (80) 

u(l,  t)  =U  cos(cot),  (81) 


where  U  is  the  amplitude  of  the  oscillation  and  co  is  the  frequency  of 
the  oscillation  on  the  right  boundary.  The  exact  solution  is  given  by 


uexac‘(x,  t)  =  Real 


_  pA 2X 

——Ue 
—  eA2 


'll, 2 


a±  \/a2  +  4icov 
2v 


(82) 


where  i  =  \/2T . 

We  solved  the  above  two  problems  with  the  first-order  hyper¬ 
bolic  advection-diffusions  equation  given  as  in  Eq.  (2)  and  (3). 
For  each  physical  time,  we  reduced  the  residuals  by  two  orders 
of  magnitude  before  advancing  in  time.  During  each  time  step, 
we  also  relaxed  the  linear  system  using  GS  relaxations  until  two 


Third-Order  RD  Schemes 


Fourth-Order  RD  Schemes 


Sixth-Order  RD  Schemes 


(a)  third-order  (Re  =  1) 


(b)  fourth-order  (Re  =  1) 


(c)  sixth-order  (Re  =  10) 


Fig.  8.  L,  error  of  the  proposed  high-order  RD  hyperbolic  schemes  for  the  boundary  layer  problem  on  nonuniform  grids. 
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Table  7 

Spatial  accuracy  with  the  third-order  RD-GT  scheme  for  the  boundary  layer  problem 
with  Re  =  1 .0. 


Number  of  nodes 

Li  error  of  u 

Order 

L\  error  of  p 

Order 

10 

9.65E-03 

8.68E-02 

25 

4.70E-04 

3.30 

4.47E-03 

3.24 

50 

4.94E-05 

3.25 

4.70E-04 

3.25 

100 

5.56E-06 

3.15 

5.26E-05 

3.16 

200 

6.58E-07 

3.08 

6.16E-06 

3.09 

300 

1.9  IE— 07 

3.06 

1.79E-06 

3.05 

500 

4.07E-08 

3.03 

3.79E-07 

3.04 

1000 

5.04E-09 

3.01 

4.66E-08 

3.02 

2000 

6.34E-10 

2.99 

5.79E-09 

3.01 

5000 

4.89E-11 

2.80 

3.68E-10 

3.00 

Table  8 

Spatial  accuracy  with  the  fourth-order  RD-GT  schemes  for  the  boundary  layer 
problem  with  Re  =  1 .0.  Note  that  when  the  differences  between  the  numerical  results 
and  exact  values  approach  the  machine  zero,  the  /.,  slope  approaches  zeroth-order. 


Number  of  nodes 

Li  error  of  u 

Order 

Li  error  of  p 

Order 

10 

8.48E-03 

6.43E-02 

25 

2.88E-04 

3.69 

1.84E-03 

3.88 

50 

1.78E-05 

4.02 

1.04E-04 

4.14 

100 

1.06E-06 

4.07 

5.87E-06 

4.15 

200 

6.41  E-08 

4.05 

3.42E-07 

4.10 

300 

1.25E-08 

4.04 

6.57E-08 

4.07 

500 

1.61E-09 

4.01 

4.35E-09 

4.05 

1000 

1.06E-10 

3.92 

4.98E-10 

4.06 

2000 

1.42E-11 

- 

1.42E-11 

- 

5000 

9.69E-12 

- 

4.75E-11 

- 

orders  of  magnitude  reduction  in  the  linear  system  residuals  was 
achieved  (see  Figs.  9  and  10.)  We  also  note  that  more  residual 
reduction  in  pseudo  time  may  be  necessary  on  more  complex 
problems  but  two  orders  of  magnitude  reduction  in  the  residuals 
were  sufficient  for  the  problems  presented  here. 

We  examined  the  convergence  rate  of  these  problems  on  sev¬ 
eral  uniform  and  nonuniform  grid  systems.  Given  in  Tables  10 
and  1 1  are  the  average  numbers  of  GS  relaxations  per  Newton  iter¬ 
ation  obtained  over  100  time  steps  for  the  periodic  and  the  oscilla¬ 
tory  boundary  condition  problems,  respectively.  Clearly  the 
convergence  rate  of  the  GS  relaxation  is  of  O(N),  not  0(N2)  as  typ¬ 
ical  for  numerical  methods  for  the  advection-diffusion  equation. 
Observe  that  for  most  grid  systems  only  two  Newton  iterations 
were  sufficient  to  obtain  accurate  solutions  regardless  of  the 
scheme  order.  We  also  note  that  we  used  At  =  0.01  for  all  the  grid 
systems.  The  maximum  grid  spacing  used  is  about  0.04.  The  corre¬ 
sponding  CFL  value  is  then  about  6.5  (based  on  At  =  0.01 ),  which  is 
significantly  smaller  that  the  maximum-allowable  CFL  values 
obtained  with  the  Fourier  analysis  for  all  the  BDF  methods  (see 
Appendix  A).  Note  that  the  maximum-allowable  CFL  value 
increases  with  grid  refinement.  The  time  step  is  orders-of-magni- 
tude  larger  than  that  required  for  conventional  explicit  schemes, 
which  is  limited  by  0(h2).  Of  course,  conventional  implicit  schemes 
also  allow  unconditionally  large  time  steps,  but  it  requires  0(N2) 
convergence  in  an  iterative  linear  solver  and  potentially  a  much 
larger  number  of  outer  iterations  as  well  if  the  exact  linearization 
is  not  possible  and  Newton’s  method  cannot  be  constructed. 
Hence,  the  method  developed  here  has  two  major  advantages  over 
conventional  methods:  the  second  order  Jacobian  formulation  and 
O(N)  iterative  convergence  in  the  linear  solver.  The  latter  advan¬ 
tage  can  be  potentially  huge  with  increase  of  the  grid  system  as 
the  speed-up  factor  is  O(N)  and  thus  grows  for  finer  grids.  Note  also 
that  the  second-order  Jacobian  formulation  is  the  advantage  of  the 
RD  scheme  over  finite  volume  schemes  where  the  compact 
Jacobian  formulation  is  only  first-order.  The  results  also  show  that 


Table  9 

Spatial  accuracy  with  the  sixth-order  RD-GT  scheme  for  the  boundary  layer  problem 
with  Re  =  10. 


Number  of  nodes 

Lt  error  of  u 

Order 

Li  error  of  p 

Order 

10 

4.81E-02 

5.54E-01 

25 

6.89E-05 

4.75 

6.19E-04 

4.89 

50 

4.54E-07 

6.11 

3.76E-06 

6.11 

100 

3.97E-09 

6.78 

3.21E-08 

6.73 

200 

3.56E-11 

7.59 

4.28E-10 

7.37 

300 

1.14E-11 

- 

1.40E-10 

- 

500 

9.23E-12 

- 

1.19E-10 

- 

the  convergence  rate  is  the  same  for  all  the  developed  high-order 
RD  schemes.  It  means  that  the  only  cost  to  these  higher  order 
schemes  is  the  evaluation  of  the  first  and  the  second  derivatives 
of  the  source  terms  and  in  the  case  of  the  fourth-order  scheme  with 
the  generalized  trapezoidal  rule,  the  second  derivatives  are  not 
required. 

We  verified  the  order  of  accuracy  of  the  proposed  schemes  on 
time-dependent  linear  problems  with  consistent  space-time  dis¬ 
cretization;  i.e„  third-order  with  BDF3,  fourth-order  with  BDF4, 
and  sixth-order  with  BDF6.  The  same  spatial  order  of  accuracy 
can  be  observed  with  the  A-Stable  BDF2  with  small  enough  time 
steps  such  that  the  temporal  error  is  comparable  to  the  spatial 
error  (i.e.,  At2  ~  hm,  where  m  is  the  scheme  order).  But  the 
consistent  pair  of  BDF  and  spatial  discretization  allows  about  two 
orders  of  magnitude  reduction  in  the  number  of  time  steps  for 
the  finest  grid  cases.  Note  that  time-accurate  computations  are 
started  by  BDF1  in  the  first  time  step  with  extremely  small  time 
step  (e.g.  At  =  10  s),  and  then  by  higher  order  BDFs  thereafter  with 
much  larger  time  steps  (e.g.  At  =  0.01).  This  will  ensure  the  order 
of  accuracy  of  the  developed  higher  order  schemes  through  all 
times.  We  remark  that  explicit  time  stepping  is  not  available  for 
time-accurate  computations  with  the  hyperbolic  system  method 
(see  Ref.  [6]  for  more  details). 

Fig.  11  shows  the  L i  error  convergence  for  the  above  time- 
dependent  linear  problem  with  the  oscillatory  boundary  condition 
(see  also  Tables  12-14),  where  clearly  the  order  of  accuracy  of  the 
proposed  schemes  are  verified  for  the  linear  time-dependent 
problems.  The  results  were  obtained  at  t  =  1.0.  We  used  the  same 
At  among  the  third-,  fourth-  and  sixth-order  schemes  and  were 
able  to  obtain  the  desired  order  of  accuracy.  (Note  that  we  showed 
the  RD  schemes  that  were  developed  with  the  generalized 
trapezoidal  rule  as  these  are  more  efficient  and  less  expensive  than 
the  ones  developed  with  divergence  formulation  of  the  source 
terms.) 

6.4.  Unsteady  nonlinear  advection-diffusion 

Consider  the  unsteady  nonlinear  viscous  Burgers  equation  with 
an  unsteady  time-dependent  source  term: 

d,u  +  dj  =  dx{vux)+S(x,t),  X  e(0,l),  (83) 

where  /  =  u2/ 2,  v  =  u,  and 

S(x,  t)  =  dtue +  ^dx(ue)2  -  d x(uepe),  (84) 

where  pc  =  dxue.  The  source  term  has  been  generated  by  the  follow¬ 
ing  function: 

(sinh  (xJico/v)  \ 

- v.  - j2L  uelmt  I  -|-  c,  C>1.  (85) 

sinh  ( s/ico/v)  ) 

so  that  it  is  the  exact  solution  to  Eq.  (83)  with  the  boundary 
conditions  defined  as 
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(a)  u(x,t  =  1.0)  (b)  p(x,t  =  1.0)  =  du/dx 

Fig.  9.  Time-dependent  linear  advection-diffusion  problem  (Re  =  33.33)  with  periodic  boundary  condition  on  N  =  10  uniform  nodes  (At  =  0.01.) 


x  x 

(a)  u(x,t  =  1.0)  (b)  p(x,t  =  1.0)  =  du/dx 

Fig.  10.  Time-dependent  linear  advection-diffusion  problem  with  oscillatory  boundary  condition  (co  =  7n/2,  k  =  2n,  U  =  1,  v  =  l)onN  =  10  uniform  nodes  (At  =  0.01.) 


Table  10 

Unsteady  linear  advection-diffusion  problem  with  periodic  BC  (a  =  l,v  =  0.03)  on 
uniform  grids.  Average  data  over  100  time  steps  are  given  (convergence  criteria: 
residuals  <  10-2). 


Number  of 
nodes 

RD-GT  scheme 
order 

GS  relaxations/Newton 
iteration 

Newton 

iteration 

25 

3rd 

5 

4 

4th 

5 

4 

6th 

6 

5 

100 

3rd 

24 

2 

4th 

24 

2 

6th 

24 

2 

300 

3rd 

33 

2 

4th 

35 

2 

6th 

35 

2 

500 

3rd 

55 

2 

4th 

55 

2 

6th 

55 

2 

1000 

3rd 

116 

2 

4th 

116 

2 

6th 

116 

2 

Table  11 

Unsteady  linear  advection-diffusion  problem  with  oscillatory  BC  (co  =  7n/2,a  =  1.) 
on  nonuniform  grids.  Average  data  over  100  time  steps  are  given  (convergence 
criteria:  residuals  <  10-2). 


Number  of 
nodes 

RD-GT  scheme 
order 

GS  relaxations/Newton 
iteration 

Newton 

iteration 

25 

3rd 

35 

4 

4th 

37 

4 

6th 

45 

4 

50 

3rd 

69 

4 

4th 

72 

4 

6th 

72 

4 

100 

3rd 

136 

4 

4th 

142 

4 

6th 

142 

4 

500 

3rd 

650 

4 

4th 

680 

4 

6th 

680 

4 

1000 

3rd 

1283 

4 

4th 

1332 

4 

6th 

1332 

4 
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Third-Order  RD  Scheme 
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Fig.  11.  Spatial  accuracy  of  the  proposed  high-< 


Fourth-Order  RD  Schemes 


(b)  fourth-order  +  BDF4 


RD-GT  hyperbolic  schemes  for  the  time-dependent 


Sixth-Order  RD  Scheme 


lo®io  h 

(c)  sixth-order  +  BDF6 


linear  problem  with  oscillatory  BC  on  uniform  grids. 


Table  12 

Spatial  accuracy  for  the  linear  time  dependent  problem  with  oscillatory  BC 
(co  =  77T/2,  a  =  1.)  using  the  third-order  RD-GT  scheme  with  the  BDF3  time 
discretization. 


Table  15 

Unsteady  nonlinear  advection-diffusion  problem  (v  =  a  =  u)  with  oscillatory  BC 
(l/  =  1,C  =  2,co  =  7n/2)  on  nonuniform  grids.  Average  data  over  1000  time  steps  are 
given  (convergence  criteria:  GS  relaxation  <  10-2 ;  Newton  residuals  <  10-8). 


Number  of  nodes 

At (BDF3) 

L]  error  of  u 

Order 

L !  error  of  p 

Order 

10 

2.50E-03 

1.40E-04 

2.86E-04 

20 

2.50E-03 

9.06E-06 

3.95 

1.90E-05 

3.91 

50 

1.25E-03 

3.54E-07 

3.54 

7.38E-07 

3.54 

100 

5.00E-04 

2.57E-08 

3.78 

4.71  E-08 

3.97 

200 

2.50E-04 

2.59E-09 

3.31 

4.40E-09 

3.42 

Table  13 

Spatial  accuracy  for  the  linear  time  dependent  problem  with  oscillatory  BC 
(co  =  7 it f2,  a  =  1.)  using  the  fourth-order  RD-GT  scheme  with  the  BDF4  time 
discretization. 

Number  of  nodes 

At (BDF4) 

Lj  error  of  u 

Order 

Lj  error  of  p 

Order 

10 

2.50E-03 

1.32E-04 

2.83E-04 

20 

2.50E-03 

7.22E-06 

4.19 

1.63E-05 

4.12 

50 

1.25E-03 

1.70E-07 

4.09 

4.01  E-07 

4.04 

100 

5.00E-04 

1.04E-08 

4.03 

2.51  E-08 

4.00 

200 

2.50E-04 

6.78E-10 

3.94 

1.64E-09 

3.94 

Table  14 

Spatial  accuracy  for  the  linear  time  dependent  problem  with  oscillatory  BC 
(co  =  771/2,  a  =  1.)  using  the  sixth-order  RD-GT  scheme  with  the  BDF6  time 
discretization. 

Number  of  nodes 

At (BDF6) 

Z-!  error  of  u 

Order 

Li  error  of  p 

Order 

10 

2.50E-03 

3.77E-06 

2.29E-05 

20 

2.50E-03 

7.19E-08 

5.71 

3.27E-07 

6.13 

50 

1.25E-03 

5.81E-09 

6.20 

3.70E-08 

5.37 

100 

5.00E-04 

3.62E-10 

5.43 

1.34E-09 

6.50 

u(0,t)=C,  (86) 

u(l,t)=C+Ucos(cot),  (87) 


where  a>  is  the  frequency  of  the  oscillation  on  the  right  boundary, 
and  U  is  the  amplitude  of  the  oscillation.  We  note  that  the  constant 
C  must  be  greater  than  1  in  order  for  the  diffusion  coefficient  to  be 
positive.  We  solved  this  time-dependent  nonlinear  advection- 
diffusion  equation  with  the  following  equivalent  first-order 
hyperbolic  system  (see  Ref.  [6]  for  more  details): 

dxu  +  dx(u2)  =  dxp  -  dcu  +  S(x,  t),  (88) 

■jdTp  =  (dxti-p/v).  (89) 


Number  of 
nodes 

RD-GT  scheme 
order 

GS  relaxations/Newton 
iteration 

Newton 

iteration 

50 

3rd 

435 

10 

4th 

430 

10 

6th 

431 

10 

100 

3rd 

879 

10 

4th 

868 

10 

6th 

864 

10 

200 

3rd 

1772 

10 

4th 

1749 

10 

6th 

1737 

10 

For  this  nonlinear  unsteady  problem,  the  manufactured  source  term 
contains  terms  that  are  already  in  the  divergence  form;  i.e.  the 
residual  evaluation  of  these  terms  are  exact.  The  d,ue  term  is  the 
only  term  in  the  manufactured  source  term  with  a  non-exact  resid¬ 
ual  evaluation.  The  BDF  discretization  of  the  dtu  term  in  Eq.  (88)  will 
not  be  in  the  divergence  form  and  therefore  will  not  have  an  exact 
residual  evaluation.  In  addition,  the  second  equation  has  a 
nonlinear  source  term,  p/v  (note  that  here  v  =  it).  We  obtained 
the  high-order  results  by  evaluating  all  of  these  non-exact  residuals 
using  the  proposed  techniques.  Newton  iterations  are  taken  to  be 
converged  when  the  overall  residuals  are  dropped  by  eight  orders 
of  magnitude.  For  each  Newton  iteration,  we  relaxed  the  linear  sys¬ 
tem  until  the  residuals  are  reduced  by  two  orders  of  magnitude. 

The  O(N)  convergence  rate  of  the  GS  relaxations  was  once  again 
achieved  for  the  time-dependent  nonlinear  hyperbolic  advection- 
diffusion  system.  This  is  given  in  Table  15,  where  the  average 
number  of  iterations  were  obtained  over  1000  time  steps  (over 
17  periods).  Note  also  that  the  high-order  RD  schemes  converged 
with  only  ten  Newton  iterations  with  the  compact  second-order 
Jacobian  formulation.  It  is  remarkable  that  such  a  rapid  conver¬ 
gence  is  achieved  for  all  high-order  schemes  with  the  second-order 
Jacobian. 

Shown  in  Fig.  12  are  the  order  of  accuracy  plots  obtained  for  this 
unsteady  advection-diffusion  problem  on  series  of  nonuniform 
grids  (see  also  Tables  16-18).  The  results  confirm  the  high-order 
accuracy  of  developed  RD  schemes  for  unsteady  nonlinear  prob¬ 
lems  on  nonuniform  grids.  We  also  observe  that  our  proposed 
sixth-order  RD  scheme,  as  discussed  in  the  previous  section,  have 
in  fact  produced  seventh-order  accurate  results  demonstrating  its 
power  characteristics  in  efficiently  providing  very  high  accurate 
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Third-Order  RD  Scheme 


(a)  third-order  +  BDF3 


Fourth-Order  RD  Schemes 


Sixth-Order  RD  Scheme 
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Fig.  12.  Spatial  accuracy  of  the  proposed  high-order  RD-GT  hyperbolic  schemes  for  the  time-dependent  nonlinear  problem  (v  =  a  =  u)  with  oscillatory  BC 
(U  =  1,  C  =  2,  to  =  77t/2)  on  nonuniform  grids. 


Table  16 

Spatial  accuracy  for  the  time-dependent  nonlinear  problem  (v  =  a  =  u)  with  oscilla¬ 
tory  BC  (11  =  1,  C  =  2,  to  =  7n/2)  using  the  third-order  RD-GT  scheme  and  BDF3 
time  discretization  on  nonuniform  grids. 


Number  of  nodes 

At  (BDF3) 

L]  error  of  u 

Order 

L\  error  of  p 

Order 

25 

1.25E-03 

6.83E-05 

3.59E-04 

50 

1.25E-03 

3.68E-06 

4.21 

1.79E-05 

4.33 

100 

1.00E-03 

2.33E-07 

3.98 

1.12E-06 

4.00 

200 

5.00E-04 

1.70E-08 

3.78 

8.23E-08 

3.77 

500 

2.50E-04 

9.58E-10 

3.14 

5.03E-09 

3.05 

Table  17 

Spatial  accuracy  for  the  time-dependent  nonlinear  problem  (v  =  a  =  u)  with  oscilla¬ 
tory  BC  (U  =  1 .  C  =  2,  w  -  771/2)  using  the  fourth-order  RD  schemes  and  BDF4  time 
discretization  on  nonuniform  grids. 


Number  of  nodes 

At  (BDF4) 

L !  error  of  u 

Order 

Lt  error  of  p 

Order 

25 

1.25E-03 

6.75E-05 

3.49E-04 

50 

1.25E-03 

6.21E-07 

4.27 

2.71E-06 

4.23 

100 

1.00E-03 

1.85E-07 

4.24 

7.68E-07 

4.40 

200 

1.00E-03 

1.07E-08 

4.11 

3.69E-08 

4.38 

500 

5.00E-04 

3.14E-10 

3.85 

6.33E-10 

4.44 

Table  18 

Spatial  accuracy  for  the  time-dependent  nonlinear  problem  (v  =  a  =  u)  with  oscilla¬ 
tory  BC  (LI  =  1,  C  =  2,  co  =  771/2)  using  the  sixth-order  RD-GT  scheme  and  BDF6 
time  discretization  on  nonuniform  grids. 


Number  of  nodes 

At  (BDF6) 

Li  error  of  u 

Order 

Li  error  of  p 

Order 

50 

1.25E-03 

7.99E-08 

5.61  E-08 

75 

5.00E-04 

5.42E-09 

6.60 

3.69E-09 

6.71 

100 

5.00E-04 

7.54E-10 

6.85 

5.09E-09 

6.91 

200 

5.00E-04 

4.30E-11 

7.06 

2.83E-09 

7.11 

solutions  and  gradients  (see  Table  18  for  more  details).  We  note, 
however,  that  the  seventh-order  accuracy  is  not  a  general 
feature  of  the  sixth-order  scheme.  Although  it  is  not  clear  from 
the  exact  solution,  it  could  be  due  to  the  sixth-order  derivative  of 
the  source  term  in  the  leading  truncation  error  being  vanishingly 
small.  We  also  note  that  the  CFL,  corresponding  to  the  largest  grid 
spacing  of  about  0.04  used  in  this  study,  is  about  30  for  the 
At  =  0.00125.  This  value  is  still  within  the  range  of  the  stability 
of  all  the  BDF  methods  (see  Appendix  A  for  more  details.)  Note 
again  that  the  maximum-allowable  CFL  increases  with  grid 
refinement. 


7.  Conclusions 

In  this  paper,  we  have  developed  a  series  of  efficient  high-order 
Residual-Distribution  (RD)  schemes  for  general  advection-diffu- 
sion  problems.  Third-  and  fourth-order  RD  schemes  were 
developed  with  the  divergence  formulation  of  the  source  term. 
These  schemes  are  very  economical  as  they  require  only  the 
evaluation  of  the  first  and  second  derivatives  of  the  source  term, 
and  not  of  the  solution.  The  first  and  second  derivatives  need  to 
be  second  and  first  order  accurate,  respectively,  on  arbitrary  grids, 
and  they  are  obtained  by  a  compact  quadratic  fit.  Third-,  fourth-, 
and  sixth-order  RD  schemes  are  also  developed  with  a  generalized 
trapezoidal  rule.  The  third-order  schemes  require  the  evaluation  of 
the  first  and  second  derivatives  of  the  source  term,  which  must  be 
second  and  first  order  accurate,  respectively.  On  the  other  hand, 
the  fourth-order  scheme  requires  only  the  second-order  accurate 
gradients,  and  therefore  is  the  least  expensive  scheme  among  all 
the  developed  high-order  schemes  in  this  paper.  All  third-  and 
fourth-order  schemes  are  constructed  within  a  five-point  stencil 
in  the  interior,  a  four-point  stencil  at  the  nodes  adjacent  to  the 
boundary,  and  a  three-point  stencil  at  the  boundary  nodes.  For 
the  sixth-order  scheme,  the  evaluation  of  the  first  and  second 
derivatives  of  the  source  term  is  required,  and  they  must  be  third- 
and  second-order  accurate,  which  is  achieved  by  a  cubic  fit.  It 
results  in  the  stencil  extension  by  one  or  two  nodes  in  the  nodal 
residual.  In  addition,  the  analysis  of  the  proposed  sixth-order  RD 
scheme  as  well  as  the  results  presented  here  show  that  this 
high-order  RD  scheme  can,  in  some  cases,  produce  seventh-order 
results  on  nonuniform  grids.  An  implicit  steady  solver  is 
constructed  based  on  the  Jacobian  derived  from  the  compact  sec¬ 
ond-order  RD  scheme.  It  is  demonstrated  that  the  solver  achieves 
a  rapid  convergence  like  Newton’s  method  for  all  high-order 
schemes  despite  the  fact  that  the  Jacobian  is  not  exact.  Specifically, 
it  requires  only  a  small  number  of  Newton  iterations,  typically  less 
than  10,  for  both  steady  and  unsteady  problems,  even  for  highly 
refined  grids,  up  to  100,000  nodes.  We  have  demonstrated  also 
that  all  of  these  high-order  schemes  are  capable  of  producing  both 
high-order  accurate  solution  and  gradient  on  nonuniform  grids 
very  efficiently  by  less  than  ten  Newton  iterations. 

The  study  presented  in  this  paper  should  be  of  interest  to 
researchers  working  on  finite-volume  schemes  because  the  RD 
scheme  is  known  to  be  equivalent  to  the  finite-volume  scheme  in 
one  dimension  (see  for  example  Ref.  [24]).  Specifically,  a  second- 
order  upwind  RD  scheme  is  equivalent  to  the  first-order  finite-vol¬ 
ume  scheme  with  a  special  form  of  source  term  discretization  [2], 
It  implies  that  the  developed  high-order  RD  schemes  may  also  be 
implemented  in  the  form  of  first-order  finite-volume  schemes  with 
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special  source  term  discretization  formulas.  The  resulting  finite- 
volume  schemes  will  be  different  from  many  other  finite-volume 
schemes  in  that  they  do  not  require  computations  of  solution 
gradients. 

The  developed  high-order  schemes  could  well  bring  significant 
improvements  to  the  numerical  methods  for  practical  problems 
such  as  material  thermal  response  calculations  of  thermal 
protection  systems  of  atmospheric  entry  vehicles  [17-19],  and 
the  experimental  aeroheating  data  reduction  [20,21],  which  are 
based  on  one-dimensional  analyses  and  routinely  used  in 
industries  (e.g.  NASA).  A  particularly  useful  scheme  would  be  the 
fourth-order  scheme  based  on  the  generalized  trapezoidal  rule 
(RD-GT)  because  it  requires  only  the  second-order  accurate 
gradients  of  the  source  term.  Application  to  these  practical 
problems  should  be  undertaken  and  is  left  as  future  work. 

Extensions  to  higher  dimensions  are  highly  desired.  To  extend 
the  developed  high-order  schemes  to  higher  dimensions,  it  is 
necessary  to  employ  a  high-order  quadrature  formula  for  integrat¬ 
ing  the  flux  divergence  term,  which  has  been  integrated  exactly  in 
one  dimension  but  cannot  be  integrated  exactly  in  higher 
dimensions.  For  the  source  term  discretization,  the  divergence 
formulation  can  be  extended  relatively  straightforwardly  while  a 
discretization  formula  such  as  the  generalized  trapezoidal  rule 
remains  to  be  found.  In  particular,  a  third-order  version  is  expected 
to  be  practical  in  multi-dimensions:  a  high-order  quadrature  for¬ 
mula  with  added  edge-midpoints  along  with  a  quadratic  fit  for 
first-derivatives.  The  solution  at  the  midpoint  can  be  obtained  by 
the  Hermite  interpolation  along  each  edge,  and  thus  does  not  need 
to  be  stored  [12],  As  a  result,  the  number  of  degrees  of  freedom 
remains  the  same  as  the  second-order  scheme,  it  is  expected  to 
be  a  very  efficient  scheme  compared  with  other  high-order  meth¬ 
ods  such  as  Discontinuous  Galerkin  or  Spectral  Volume  methods. 

Finally,  extensions  to  more  complex  nonlinear  equations  such 
as  the  Navier- Stokes  equations  remain  as  a  challenge.  For  the  com¬ 
pressible  Navier- Stokes  equations,  the  complete  eigen-structure  of 
the  whole  system  has  yet  to  be  found.  The  construction  of  the 
upwind  RD  scheme  based  on  a  single  hyperbolic  system,  as  pre¬ 
sented  in  this  paper,  is  a  challenge.  To  overcome  the  difficulty,  a 
simplified  approach  has  been  proposed  and  demonstrated  for  a 
finite-volume  scheme  in  Ref.  [3  ,  which  is  based  on  the  indepen¬ 
dent  treatment  of  the  inviscid  and  viscous  terms.  A  similar 
approach  may  become  necessary  for  the  RD  schemes. 


The  spatial  operator  is  defined  as: 

^  +  B  Ja>(() .  (A.3) 

where 


-a(eip  -  1)  v(e,7l-l) 

(e,7,-l  )/Tr  %  J’ 

-a(l  -  e~'P)  v(l  -  e~'p) 
(1  -  e-'^/Tr  S^R 


(A.4) 


The  Sp  is  associated  with  the  source  term  and  therefore  depends  on 
the  proposed  source  term  discretization.  For  discussion  purposes, 
here  we  only  consider  high  order  RD-GT  schemes  and  the 
corresponding  source  terms  become: 


(  eip  +  1  +  Ar  +  Br  :  3rd-order 

s1^  =  w  <  &  +  1  +  A*  :  4th-order  (A.  5) 

(  e'P  +  1  +  CR  +  Dr  :  6th-order 

(  e~ip  +  1  +  Al  +  Bl  :  3rd-order 
=  ZT  |  e  +  1  +aL  '■  4th-order  (A.6) 

[  +  1  +  CL  +  Dl  :  6th-order 


where, 

AR  =  (e'7*  -  e“l7i  -  e2'7*  + 1)/12  (A.7) 

Br  =  (-e,7!  +  e-ip  -  e2ip  -  3)/1000  (A.8) 

CR  =  (~e~2ip  +  3e^  -  2eip  +  3e2ip  -  e3'"  -  2)/12  (A.9) 

Dr  =  (~eip  +  e~ip  +  e2ip  -  l)/60  (A.10) 

Al  =  (~eip  +  e-ip  -  e-2it>  -  1)/10  (A.ll) 

B1  =  (— e'7*  -  3e-,7J  +  e  2i/l  -  3)/ 1000  (A.  12) 

CL  =  (4e-2l7i  -  2 e-V  +  3eip  -  e2'7*  -  -  2)/12  (A.13) 

D1  =  (e1^  -  e~®  +  e~2f/i  + 1 ) /60  (A.14) 

The  mass  matrix  N/f  is  defined  as: 

N/'  =  J(B  J^.  +  B  J^),  (A.15) 

where 
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Appendix  A.  BDF  stability  analysis  on  the  proposed  high  order 
RD  schemes 


k  =  Tr 


cl> 


0 


0 

0 


(A- 16) 


Note  that  only  the  variable  u  evolves  with  the  BDF.  Thus,  we  can 
eliminate  the  p0  by  solving  the  second  equation  in  the  system 
(A.2)  and  substituting  it  back  to  the  first  equation  and  arrive  at: 

BDF(tio)  = /lio,  (A.  17) 


Consider  a  Fourier  mode  of  non  dimensional  wave  number 
P  6  [0, 7t]  on  a  uniform  mesh  of  spacing  h: 

U^  =  U0  (t)el7,<*-*f>7\  (A.l) 

where  U/f  =  (u^p11)  and  U0(t)  =  (u0(t),p0(t)).  Applying  this  Fourier 
mode  to  the  discretized  first  order  hyperbolic  system,  we  arrive  at: 

(A.  2) 

where  M1’  and  N/(  are,  respectively,  the  spatial  operators  for  the 
spatial  and  temporal  terms  of  the  system.  Note  that  the  spatial 
operator  arises  from  the  temporal  term  because  the  physical  time 
derivative  is  discretized  in  space  and  distributed  to  the  nodes  in 
our  schemes.  The  matrix  N,f  may  be  called  the  mass  matrix. 


where  BDF(u0)  denotes  the  time  discretization  of  du0/dt  by  the  BDF, 
and 


/  = 


K,i> 


Np 


(A- 18) 


With  the  physical  time  step  defined  as  At  =  CFL/i/(a  +  v/h),  the 
quantity  relevant  to  the  stability  2/*At  can  be  shown  to  depend  only 
on  two  parameters:  ReLr  =  aLr/v  and  the  mesh-Reynolds-number 
Reh  =  ah/v.  The  former  is  essentially  equivalent  to  the  global 
Reynolds  number,  Re  =  a/v  since  Lr  =  0(1)  (see  [6]  for  details);  it 
plays  a  role  of  properly  weighting  the  upwind  advection  and  the 
upwind  diffusion  schemes  [2], 

We  now  numerically  perform  stability  analysis  for  BDF3, 
BDF4,  and  BDF6  and  estimate  the  maximum-allowable  CFL 
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Re(X.) 

(a)  BDF3  with  3rd-order  RD-GT 


Re(X-) 

(b)  BDF4  with  4th-order  RD-GT 


Re(k) 

(c)  BDF6  with  6th-order  RD-GT 


Fig.  A.13.  Maximum-allowable  CFL  for  the  proposed  high-order  time-dependent  RD  schemes  (h  =  0.1 .  Re  =  100).  The  stable  region  is  bounded  by  a  curve  represented  by  the 
red  circles.  The  eigenvalues  of  the  discretizaiton  are  plotted  in  blue  and  all  contained  in  the  stable  region  for  the  maximum  CFL  number. 


Table  A.19 

Maximum-allowable  CFL  values  for  the  proposed  high-order  unsteady  RD  hyperbolic  advection-diffusion  schemes. 


Grid  spacing 

Re  =  1 

Re  =  100 

Re  =  106 

BDF3 

BDF4 

BDF6 

BDF3 

BDF4 

BDF6 

BDF3 

BDF4 

BDF6 

0.1 

1000 

500 

20 

1.5 

0.5 

0.02 

0.01 

90,000 

50,000 

1800 

25 

10 

0.35 

0.001 

1500 

500 

20 

0.04 

0.005 

0.0002 

0.0001 

150,000 

50,000 

2000 

0.15 

0.05 

0.002 

0.00001 

1.5 

0.5 

0.02 

0.000001 

25 

10 

0.3 

(=  A t(a  +  v/h)/h)  values  for  different  mesh  sizes  and  Re  —  a/v 
numbers.  Shown  in  Fig.  A.13  are  the  stability  regions  and 
maximum  CFL  values  graphically  illustrated  corresponding  to  the 
proposed  time-dependent  high-order  RD-GT  schemes.  Table  A.19 
provides  maximum  CFL  for  ranges  of  Re  and  grid  spacing  for  the 
high  order  schemes.  Observe  that  the  maximum-allowable  CFL 
number  varies  significantly  with  Re  on  a  given  mesh,  but  it  is 
essentially  the  same  for  the  same  Reh.  Note  that  Reh  <  2  is  required 
to  avoid  numerical  oscillations  [2],  and  the  maximum-allowable 
CFL  number  in  this  range  is  large  enough  to  perform  practical 
computations. 
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